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}]