mathematica 10 not showing numerical solution of differential equations?

I just got the new mathematica version 10 and tried to solve the following system of differential equation.

$$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$$ $$f(0)=h(0)=0, \quad f(\infty)=h(\infty)=1$$

and my code in mathematica,

Clear[eqn, bc, f, h]
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};
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}]


Mathematica returned the error message:

NDSolve::ndsz: At r == 16.162462464292577, step size is effectively zero; singularity or stiff system suspected. >>


Now, my problem is the New version 10 does not give me the InterpolatingFunction

while the previous version 8 gives me,

{{h->InterpolatingFunction[{{1.*^-6,20.}},"<>"],f->...}}


Does anybody know what's going wrong with the new version? How may I get the solution in V10? My sincere thanks for you patience and help!

It seems that the numerical solution behaves strangely in large $r$ region. If I set the boundary at $r=50$, the solution looks like this,

• can you show a plot of the solution you obtained for f(r) and h(r) in V8 for r=0..20? What options did you use with NDSolve for V8? Same as shown? or did you change something? Jan 12 '15 at 10:40
• I get the warning in V9.0.1, but it still produces a solution: !Mathematica graphics Jan 12 '15 at 11:17
• Thank you @Michael E2 , that's exactly what I obtained from V8.0. But V10 gave me nothing Jan 12 '15 at 11:56
• Maple solves this, using BVP method with mid-point, used it since there is singularity at t=0. Here is screen shot. !Mathematica graphics Tried all the mid points methods in M with NDSolve, but none of them worked. same error. V 10.02 Jan 12 '15 at 12:26
• The output I got from V8 is the interpolatingFunction which I can use them {{h->InterpolatingFunction[{{1.*^-6,20.}},"<>"],f->...}} when I plot them out I just got the same plot as @Michael E2. However V10 does not give me that Jan 12 '15 at 12:29

The problem is that NDSolve is using the shooting method and for some of the initial conditions it tests lead to a singularity. Indeed there seems to be a small neighborhood around the desired solution where the integration does not blow up before r == 20. I don't know why previous versions of Mathematica were able to push through and find an adequate solution. In any case V10 seems to give up when it gets an error.

One approach is to transform the differential equation to be a bit better behaved near r == 0. We can replace f[r] by r^2 g[r] and h[r] by r j[r].

Clear[eqn, bc, f, h, g, j]
With[{r1 = 0.000001, r2 = 20},
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};
bc = {f[r1] == 0, f[r2] == 1, h[r1] == 0, h[r2] == 1};
{sol} = NDSolve[{eqn, bc} /. {f -> (#^2 g[#] &),
h -> (# j[#] &)}, {g, j}, {r, r1, r2}]
]

Plot[{r j[r], r^2 g[r]} /. sol // Evaluate, {r, 1.*^-6, 20}]
`

• In fact I have a second, it seems that there is some problem with the numerical solution if I set $r_2$ to a large value, say, $r_2 =50$. The numerical solution behaves strangely at large $r$ region and the shape at small $r$ also change. I will upload the screenshot in my question. Jan 13 '15 at 5:33