Solving a bad-behaving equation with shooting method

Question

I would like to numerically solve

$a^2\frac{d}{d\sigma}[(1+A_tQ)\frac{d\psi}{d\sigma}]=\psi(1+A_tQ-a\Phi\frac{dQ}{d\sigma})$

where $Q(\sigma)=\frac{2}{\sqrt{\pi}}\int_0^\sigma exp(-x^2)dx$.

$a$ and $A_t$ are certain parameters and $\psi(x)$is the desired function with certain boundary values and $\Phi$ is the eigenvalue of the equation. I started with the built-in shooting method:

a = 1/5;
A = 0.843;
Q[x_] := 2/Sqrt[Pi]*Integrate[Exp[-p^2], {p, 0, x}]; 
eq := a^2*D[(1 + A*Q[x])*D[psi[x], x], x] ==  psi[x]*(1 + A*Q[x] - a*Phi[x]*2/Sqrt[Pi]*Exp[-x^2]);
bc = {psi'[-2] == -1/a*Exp[-2/a], psi[-2] == -Exp[-2/a], psi[100] == 0};
s = NDSolve[{eq, bc, Phi'[x] == 0}, {psi[x], Phi[-2]}, {x, -2, 100}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {Phi[-2] == 3}}];

However, this equation behaves badly, and with shooting method I always get NDSolve:berr, and the solution of Phi stays where the initial condition is set.

So I wrote a simple shooting method myself:

Q[x_] := 2/Sqrt[Pi]*Integrate[Exp[-p^2], {p, 0, x}]; 
eq := a^2*D[(1 + A*Q[x])*D[psi[x], x], x] == psi[x]*(1 + A*Q[x] - a*Phi*2/Sqrt[Pi]*Exp[-x^2]);
bc := {psi'[-2] == -1/a*Exp[-2/a], psi[-2] == -Exp[-2/a]};
a = 1/5;
A = 0.843;
Phi = 3;
d = 0.1;
Sol = NDSolve[{eq, bc}, psi, {x, -2, 100}, 
   Method -> "ExplicitRungeKutta"];
m1 = Maximize[{Abs[psi[x] /. Sol[[1]]], 80 <= x <= 100}, x][[1]];
While[d > 0.0001,
 Phi = Phi + d; 
 Sol = NDSolve[{eq, bc}, psi, {x, -2, 100}, 
   Method -> "ExplicitRungeKutta"];
 m1new = Maximize[{Abs[psi[x] /. Sol[[1]]], 80 <= x <= 100}, x][[1]];
 If[m1new > m1, Phi = Phi - d; d = d/2, m1 = m1new];
 Print[Phi];
 ]

This just works better and gave an approximate value, although the solution of $\psi(\sigma)$ is still diverging.

Since my scanning scheme is rough and doesn't guarantee a best solution in principle, is there any way to get the built-in shooting method to work? Or is there any better approach to this job?

Answer 1

To solve this problem we don't need shouting method, also range $-2\le x\le 100$ is to large since even for $-2\le x\le 5$ we should compute with WorkingPrecision -> 30. Numerical solution is given by

a = 1/5;
A = 843/1000;
Q[x_] := Erf[x];
eq = a^2*D[(1 + A*Q[x])*D[psi[x], x], x] == 
   psi[x]*(1 + A*Q[x] - a*Phi*2/Sqrt[Pi]*Exp[-x^2]);
bc = {psi'[-2] == -1/a*Exp[-2/a], psi[-2] == -Exp[-2/a]};
s = ParametricNDSolveValue[{eq, bc}, psi[5], {x, -2, 5}, {Phi}, 
   WorkingPrecision -> 30];

sol = FindRoot[s[z] == 0, {z, 378/100}]

(*Out[]= {z -> 3.78066}*)

   
Phi0 = Rationalize[z /. sol, 10^-31];
    
sol1 = NDSolveValue[{eq /. Phi -> Phi0, bc}, psi, {x, -2, 5}, 
   WorkingPrecision -> 30];

Visualization

Plot[sol1[x], {x, -2, 5}, Frame -> True, GridLines -> Automatic]

As we can see from this picture we can't extend numerical solution up to x=100 since asymptotic solution psi is proportional to $e^{-5x}$, therefore psi[100] is about Exp[-500].