# (NDSolve) Non-linear 2nd order ODE, regular singular point (looking for good methods for this problem)

I am solving this set of non-linear 2nd order ODE by NDSolve, $$r^2\frac{d^2f}{dr^2} = 2f(1-f)(1-2f)+\frac{r^2}{4}h^2(f-1)$$ $$\frac{d}{dr}\left[r^2\frac{dh}{dr} \right]=2h(1-f)^2+\lambda r^2(h^2-1)h$$ $$f(0)=h(0)=0, \quad f(\infty)=h(\infty)=1$$

From the series solution, I also know the behavior at $r\approx0$, $$f\approx \alpha r^2, \quad h\approx \beta r$$ Therefore the singularity at $r=0$ is a regular singular point.

I tried using NDSolve,

\[Lambda] = 625/2048;
eqn = {r^2*f''[r] ==
2 f[r] (1 - f[r]) (1 - 2 f[r]) - r^2/4 (h[r])^2 (1 - f[r]),
D[r^2*h'[r], r] == 2 h[r] (1 - f[r])^2 + \[Lambda]*r^2 (h[r]^2 - 1) h[r]};
bc = {f[0.000001] == 0, f[20] == 1, h[0.000001] == 0, h[20] == 1};
sol = NDSolve[{eqn, bc}, {h, f}, {r, 0.000001, 20}]


However, it doesn't solve it says At r =...., step size is effectively zero; singularity or stiff system suspected.

I also tried,

NDSolve[{eqn, bc} /. {f -> (#^2 g[#] &), h -> (# j[#] &)}, {g, j}, {r, r1, r2}]


with boundary condition at r1 and r2. Still it doesn't give me the solution.

I guess there are some subtleties about non-linear ODEs. Since it is non-linear, the solution is highly sensitive to boundary conditions. Unless I set the exact boundary condition, it won't give me the right answer. However, the ODE above has a regular singular singular point at $r=0$, which I don't know how to impose the exact boundary condition at $r=0$ in mathematica. It stops calculation immediately if I set the boundary condition at $r=0$.

I am looking for help from you guys, especially those had the same problem before. My sincere thank for your help!

Do you know if the method of relaxation help in this case?

Update: The method of relaxation works well for this type of problem, however it seems that Mathematica does not implement this for NDSolve algorithm.

• at r=0 is a regular singular point. I though the classification of ordinary, regular, and not regular singular points applies to linear ODE's. But what you have is non-linear ODE's. I could be wrong, but I remember something like this from a class. Commented Mar 23, 2015 at 4:56
• you are right but how do you describe that in nonlinear ODEs? Commented Mar 23, 2015 at 5:59
• @Nasser The classification of ordinary/regular/singular point can be applied to nonlinear ODEs by investigating the asymptotic behavior of their solution in the vicinity of that point. Commented Jun 6, 2018 at 18:01

Using the final substitution given in the Question with the addition of Method -> "StiffnessSwitching" produces a solution for r2 as large as 7.4, after which even "StiffnessSwitching" is insufficient to treat the equations at f and h very near 1.

λ = 625/2048;
r1 = 10^-6; r2 = 7.4;
eqn = {r^2*D[D[f[r], r], r] ==
2 f[r] (1 - f[r]) (1 - 2 f[r]) - r^2/4 (h[r])^2 (1 - f[r]),
D[r^2*D[h[r], r], r] ==
2 h[r] (1 - f[r])^2 + λ*r^2 (h[r]^2 - 1) h[r]};
bc = {f[r1] == 0, f[r2] == 1, h[r1] == 0, h[r2] == 1};
sol = NDSolveValue[{eqn, bc} /. {f -> (#^2 g[#] &), h -> (# j[#] &)}, {g, j},
{r, r1, r2}, Method -> "StiffnessSwitching"];
Plot[{r^2 sol[[1]][r], r sol[[2]][r]}, {r, r1, r2}]


Without "StiffnessSwitching", integration fails at r2 of about 5.22.

In the limit of large r, the first term on the right side of each equation is small compared to the other terms and can be dropped;

{r^2*D[D[f[r], r], r] == - r^2/4 (h[r])^2 (1 - f[r]),
D[r^2*D[h[r], r], r] == λ*r^2 (h[r]^2 - 1) h[r]};


Moreover, since f and h approach 1 for large r, the equations can be linearized about 1 to yield

{r^2*D[D[f[r], r], r] == r^2/4 (f[r] - 1),
D[r^2*D[h[r], r], r] == 2 λ*r^2 (h[r] - 1)};


These equations are uncoupled and can be solved analytically (or with DSolve).

fasym = 1 - a Exp[-r/2];
hasym = 1 - b Exp[-Sqrt[2 λ] r]/r


with a and b constants. (I thank mastrok for pointing out in a Comment below that I had mistakenly assumed a = b = 1 in the first version of this addendum.) These asymptotic expressions can be used to improve the boundary conditions for NDSolve. To do so, set

f[r2] - 1 == - a Exp[-r2/2];
f'[r2] == a Exp[-r2/2]/2;


a can be eliminated between these two equations to yield 2 f'[r2] + f[r2] == 1. A similar but more complicated expression is obtained for h. In all,

bc = {f[r1] == 0, 2 f'[r2] + f[r2] == 1, h[r1] == 0,
r2 h'[r2] + h[r2] (Sqrt[2 λ] r2 + 1) == 1 + Sqrt[2 λ] r2};


Although the improved boundary conditions do not increase the range over which NDSolve produces meaningful results (in fact, it reduces the upper bound from 7.44 to 6.99), it does moderately improve the overall results. The figure below shows the numerical (solid) and asymptotic (dashed) curves superimposed.

Show[Plot[{r^2 sol[[1]][r], r sol[[2]][r]}, {r, r1, r2}],
Plot[{1 + a Exp[-r/2], 1 + b Exp[-r Sqrt[2 λ]]/r}, {r, 5, 10},
PlotStyle -> Dashed], PlotRange -> All]


The excellent overlap between the small and large r solutions indicates that the asymptotic solutions are quite accurate for r as small as 5. Quantitatively, a and b are equal to 1.36 and 3.37 for r2 = 6.99 and 1.41 and 3.29 for r2 = 6.00.

• useful answer. one thing I don't get is, why are the constants in the asymptotes 1? I got 2 constants in these expressions, f = 1 - AExp[-r/2]; h = 1 - BExp[-Sqrt[2 λ] r]/r Commented Mar 24, 2015 at 4:16
• Also, the solution is extremely sentive to r1 and r2. Changing them by a bit will totally change the numerical solution. That's the annoying part for me. Is there any way to impose better boundary condition ? Commented Mar 24, 2015 at 4:31
• @mastrok I have not found the results to be sensitive to r1, and only modestly sensitive to r2 if NDSolve succeeds in satisfying the boundary conditions accurately at r2. More to follow on your first comment. Commented Mar 24, 2015 at 12:35
• Do you know if there is other better method for solving these ODEs ? Like relaxation, I just heard about this method but not sure about the details Commented Mar 26, 2015 at 10:47
• The relaxation method applied to nonlinear ODEs is described here and elsewhere on the web. To use it, you would need to discretize your equations using FEM, as described in the Mathematica FEM Users Guide. Conceptually, the relaxation method can be thought of as converting the ODEs to first order PDEs in time, the solutions of which (one hopes) relax to a steady state that is the solution to the original ODEs. Commented Mar 26, 2015 at 16:42