There are two problems with the ODE at r == 0
: (1) To solve for k''[r]
you have to divide by r^2
, that is, divide by zero, (2) and k'[0] == 0
for all solutions, which frustrates the shooting method:
ode = r^2 k''[r] + 2 r k'[r] + 2 r^2 (k'[r]^2)/(1 - k[r]) - 2 (k[r] + 1) (2 - k[r]) == 0;
Solve[ode, k''[r]]
% /. k[r] -> -1 // Apart (* initial condition k == -1 at r == 0 *)
(*
{{k''[r] ->
-2 (2 - k[r] - 2 k[r]^2 + k[r]^3 - r k'[r] + r k[r] k'[r] - r^2 k'[r]^2) /
(r^2 (-1 + k[r]))}}
{{k''[r] -> -(2 k'[r])/r - k'[r]^2}}
*)
If k'[r]
approaches a finite, non-zero value at r == 0
, then k''[r]
has a pole, which would imply k'[r]
goes to infinity; therefore k'[r]
must approach zero. As I understand how the shooting method is implemented, NDSolve
would use FindRoot
to find an initial value k'[0]
such that the solution satisfies the boundary condition k[1] == 0
. But that's going to be hard to do, if the ODE is simultaneously constraining k'[0]
to be zero. The coefficient that needs to be varied is k''[0]
. (An alternative approach, such as bbgodfrey's, is to move the starting point far enough away from r == 0
so that the numerics allow NDSolve
to hit the target condition k[1] == 0
.)
To allow NDSolve
to shoot with k''[0]
, we need a third-order equation. We can differentiate the ODE to get one, and since we know we want k[0] == -1
and k'[0] == 0
, the boundary condition k[1] == 0
gives us the third condition we need for a solution. That leaves the problem of dividing by zero.
Using the shooting method with the third-order equation solves the numerics problem, and that allows us to bring the starting value for r
close to zero.
Block[{r0 = 1*^-6},
sol1 = NDSolve[
{D[ode, r],
k[r0] == -1, k'[r0] == 3.512`16 r0, k[1] == 0},
k, {r, r0, 1}]]
In fact we can help NDSolve
out a little by using Taylor's theorem to approximate k[r0]
, k'[r0]
. Inspecting the solution above, we can see that k''[r]
is about 3.512
near r == 0
.
Block[{r0 = 1*^-6, q0 = 3.512},
sol2 = NDSolve[
{D[ode, r],
k[r0] == -1 + q0 r0^2/2, k'[r0] == q0 r0, k[1] == 0},
k, {r, r0, 1}]]
Plot[k[r] /. First@sol2, {r, 1*^-6, 1}]

If we check the boundary conditions, they look pretty good:
Block[{r0 = 1*^-6},
{{k[r0] + 1, k'[r0], k''[r0]}, {k[1]}} /. First@sol2
]
(* {{1.75604*10^-12, 3.512*10^-6, 3.51366}, {-4.08313*10^-7}} *)
NDSolve
solves a DE for the highest-order derivative. At the initial condition, it will divide by zero (i.e. byr^2
forr == 0
). $\endgroup$ – Michael E2 Jul 2 '15 at 18:41$MachineEpsilon
. $\endgroup$ – J. M.'s ennui♦ Jul 2 '15 at 21:38