NDSolve with coupled ODE's and unknown singularities

Question

I have two coupled ODEs that I am trying to solve numerically. It appears that there is a singularity in the solution to the equations which I am unsure how to get past. Both functions $\alpha$ and $\phi$ should asymptotically converge to 1 at infinity, however I can't set my boundary conditions to be more than 3.4; I get the error:

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

I have seen suggestions that I should start away from a singularity (i.e. r = 0) or change boundary conditions etc, but my boundary conditions are fixed.

How do I get around this issue?

Code below:

(* Set Constants *)
n = 1;
λ = 1;

(* Set functions *)
p[r_] = 1/r^2;
q[r_] = 1/r;

(* ODE's, which for some reason I called pde1/2 *)
pde1 = D[ϕ[r], {r, 2}] + q[r]*D[ϕ[r], {r, 1}] - p[r]*(n - α[r])^2 ϕ[r] + λ/2 (1 - ϕ[r]^2) ϕ[r] == 0
pde2 = D[α[r], {r, 2}] - q[r]*D[α[r], {r, 1}] + (n - α[r]) ϕ[r]^2 == 0

(* Min/max numbers *)
min = 0.000001;
max = 3.4; (* Define a "finite infinity" - I would like to set this higher *)

(* Boundary Conditions *)
bc = {ϕ[min] == 0, α[min] == 0, ϕ[max] == 1 - 1/max Exp[-max], α[max] == 1 - 5/max Exp[-max]};

sol = NDSolve[
       {pde1, pde2, bc}, {ϕ, α}, {r, min, max}, 
       Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}, 
       AccuracyGoal -> 15, PrecisionGoal -> 15, MaxSteps -> 20000
      ]

Plot[Evaluate[{ϕ[r], α[r]} /. sol], {r, 0, 4}]

Answer 1

In this BVP, Mathematica uses the shooting method, and you're shooting at an unstable solution. You can iteratively approach the solution, extending the interval of integration at each step, but at machine precision, you can extend it only so far.

I'll set up bc to depend on symbolic min and max, so we can fiddle with them at each iteration:

(*Boundary Conditions*)
Clear[min, max];
bc = {ϕ[min] == 0, α[min] == 0, ϕ[max] == 1 - 1/max Exp[-max],
  α[max] == 1 - 5/max Exp[-max]};

Then with max = 3.4, we get an initial solution, just as the OP does.

Block[{min = 0.000001, max = 3.4},
 sol1 = NDSolve[{pde1, pde2, bc}, {ϕ, α}, {r, min, max}]
 ]

From this point, the basic iterative step uses the previous step to determine the starting initial conditions for the shooting method. So thus:

Block[{sol0 = First@sol1, min = 0.000001, max = 7.5}, 
 sol2 = NDSolve[{pde1, pde2, bc}, {ϕ, α}, {r, min, max}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" ->
       {ϕ[min] == (ϕ[min] /. sol0), ϕ'[min] == (ϕ'[min] /. sol0),
        α[min] == (α[min] /. sol0), α'[min] == (α'[min] /. sol0)}}]
 ]

The setting max = 7.5 was determined by running NDSolve first out to 10; but it developed a singularity at 7.7, so I stepped it back a little. (Push too close to the limit and the shooting method might have problems.)

Iterate again:

Block[{sol0 = First@sol2, min = 0.000001, max = 15}, 
 sol3 = NDSolve[{pde1, pde2, bc}, {ϕ, α}, {r, min, max}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" ->
       {ϕ[min] == (ϕ[min] /. sol0), ϕ'[min] == (ϕ'[min] /. sol0),
        α[min] == (α[min] /. sol0), α'[min] == (α'[min] /. sol0)}}]
 ]

Hmm, no error but something looks wrong.

Plot[Evaluate[{ϕ[r], α[r]} /. sol3], {r, 0, 15}, GridLines -> {None, {1}}]

Well, that might be as good as you can do at a WorkingPrecision of MachinePrecision.

If you use WorkingPrecision -> 20, you get the following and probably could push the integration farther than r == 15.:

An alternative approach is to implement one's own shooting method, and adapt the initial conditions to integrate as far as possible. This was one of the methods used in my answer to Nonlinear differential equation: numerical solution.