Solving a PDE with tricky boundary conditions

I am not sure if their paper is right, after I do the change of variable, I think the equation should be \begin{align*} \frac {\partial F}{\partial x}+\left(\frac {2}{\beta}\sin^{4} \theta\right)\frac {\partial^{2} F}{\partial \theta^{2}}+\left(\left(x+\frac {2}{\beta}\sin 2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right)\frac {\partial F}{\partial \theta}=0,\quad 0<\theta<\pi,\quad x<10 \end{align*} \begin{align*} F(10,\theta)=\begin{cases} \Phi\left(\frac {10-\cot^{2}\theta}{\sqrt{\left(4/\beta\right)\cot \theta}}\right)\quad &0\leq \theta\leq \pi/2\\ 1\quad &\pi/2\leq \theta\leq \pi \end{cases} \end{align*} \begin{equation*} F(x,0)=0,\quad x<10 \end{equation*} where we use $-\cot(\pi)=\infty$.

However, they might do the extension to $F(x,\theta)$ and force the $F(x,\theta)$ is a constant in $\theta$ when $\theta>\pi$. Anyway, if you want to numerically solve the following equation: \begin{align*} \frac {\partial F}{\partial x}+\left(\frac {2}{\beta}\sin^{4} \theta\right)\frac {\partial^{2} F}{\partial \theta^{2}}+\left(\left(x+\frac {2}{\beta}\sin 2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right)\frac {\partial F}{\partial \theta}=0,\quad 0<\theta<\infty,\quad x<10 \end{align*} \begin{align*} F(10,\theta)=\begin{cases} \Phi\left(\frac {10-\cot^{2}\theta}{\sqrt{\left(4/\beta\right)\cot \theta}}\right)\quad &0\leq \theta\leq \pi/2\\ 1\quad &\pi/2\leq \theta \end{cases} \end{align*} \begin{equation*} F(x,0)=0,\quad x<10 \end{equation*}

You can use the finite-difference method.

To use the finite-difference method, we first need to set up grid points: For example, you can define $$ x^0=10,\quad x^n=x^0-n\Delta x=10-n\Delta x $$ and $$ \theta_0=0,\quad \theta_m=\theta_m+m\Delta\theta\,. $$ where $m,n\in\mathbb{N}$. Then we denote the approximation of $F(x^n,\theta_m)$ by $F^n_m$.

Since you want solution at $x=-10$ with $\theta\in[0,2\pi]$, we denote denote $M=\frac{2\pi}{\Delta \theta}$ the number of grid points and $N=\frac{20}{\Delta x}$ the steps you need to run. Because it only has one side boundary condition $\theta=0$, we should use forward difference. More specifically, from step $n$ to $n+1$, according to the equation, we can get $F^{n+1}_m$ by solving \begin{equation} \begin{aligned} &\frac{F^{n+1}_m-F^{n}_m}{\Delta x}=\left(\frac {2}{\beta}\sin^{4} \theta_m\right)\frac {F^n_{m+1}-2F^n_m+F^n_{m-1}}{\Delta \theta}\\ &+\left(\left(x^n+\frac {2}{\beta}\sin 2\theta_m\right)\sin^{2}\theta_m-\cos^{2}\theta_m\right)\frac {F^n_{m+1}-F^n_m}{\Delta \theta}=0 \end{aligned}\,,\quad 0<m<M,\ 0\leq n\leq N-1\tag{1} \end{equation} And we always let $F^{n+1}_0=0$. However, there is a problem. From (1), you can notice that we haven't updated $F^{n+1}_M$. For example, you don't know how to calculate $F^1_{M}$ since you need to know $F^0_{M+1}$ if you want to use (1).

To overcome this problem, we need to use the initial condition is defined for all $\theta\geq 0$. According to the initial condition, we know the value of $$ F^0_m,\quad \forall 0\leq m\leq N+M\,. $$ Then, in each step, instead of (1), we can use \begin{equation} \begin{aligned} &\frac{F^{n+1}_m-F^{n}_m}{\Delta x}=\left(\frac {2}{\beta}\sin^{4} \theta_m\right)\frac {F^n_{m+1}-2F^n_m+F^n_{m-1}}{\Delta \theta}\\ &+\left(\left(x^n+\frac {2}{\beta}\sin 2\theta_m\right)\sin^{2}\theta_m-\cos^{2}\theta_m\right)\frac {F^n_{m+1}-F^n_m}{\Delta \theta}=0 \end{aligned}\,,\quad 0<m<M+N-n,0\leq n\leq N-1\tag{2} \end{equation} where we notice that $m$'s upper bound is changed to $M+N-n$. The reason is similar to before. Since we know $F^0_m$ for $0\leq m\leq M+N$. Using (1), we can obtain $F^1_m$ for $0\leq m\leq M+N-1$. Iteratively, we can obtain $F^n_m$ for $0\leq m\leq M+N-n$. Using (2), we can always update $F^{n}_M$ when $n\leq N$. Thus, we can obtain an approximation to $F(-10,\theta)$ for $\theta\in[0,2\pi]$. I think their paper also means similar idea by focusing on a larger interval.