I am interested in numerically solving the ODE below for the function $\theta(r)$:
$$-\theta''=\frac{1}{r}\theta'-\frac{1}{r^2}\sin(\theta)\cos(\theta)-\frac{D}{Ar}\sin^2(\theta)+\frac{K}{A}\sin(\theta)\cos(\theta)-\frac{HM_s}{2A}\sin(\theta)$$
with the boundary conditions
$$\theta(0)= 0$$ $$\theta'(r=R) = \alpha$$
where $R$ is some constant.
I have used numerical methods such as the Runge Kutta method and general ODE solvers (such as NDSolve from mathematica) before. I am not asking for help on general numerical methods.
This problem is causing me issues because at $r=0$, several terms become undefined because of the dependence on $1/r$ or $1/r^2$. But also, $\theta \rightarrow 0$, so $\sin(\theta)\rightarrow 0$, so it is possible that the orders of the $0$'s "cancel out" in terms like $-\frac{1}{r^2}\sin(\theta)\cos(\theta)$. For physical reasons, I expect that $\theta''(0)$ is finite.
So for the first timestep starting at $r=0$, the solver cannot determine $\theta''$. When the ODE is put into a solver such as Mathematica's NDSolve, the error
Power: Infinite expression 1/0.^2 encountered.
Are there any techniques for numerically solving ODE's with $1/r^n$ terms when $0$ is in the domain of interest for the ODE?
This equation comes from the paper DOI: 10.1103/PhysRevB.80.134410, where the solution for $\theta$ was numerically determined on the domain $[0, R]$. However, the paper is not clear on the details of the solution, only saying a Runge Kutta algorithm was used to solve the ODE.
Thank you for any help you can offer.
You can replace the initial segment with some other approximation, for instance a Frobenius series. For small $r$ the equation is approximately $$ r^2θ''+rθ'-θ=0, $$ using the small-angle approximation. This Euler-Cauchy DE has basis solutions $r$ and $r^{-1}$, only the first is compatible with the boundary conditions.
Insert $θ=c_1r$ or add some higher-degree terms $θ=c_1r+c_3r^3$ and evaluate the Taylor expansions to compute $c_3$ from $c_1$. Taking the linear term and using RK4 with step size $h$ after the first segment, the order $h^4$ can be obtained globally if the first segment is of length $h^2$ (probably $h^{1.5}$ is sufficient).
For the boundary value problem on $[A,R]$ with $A>0$ use $Aθ'(A)-θ(A)=0$ as left boundary condition, continue left with $c_1=θ'(A)$.