Solving Integro-differential equation numerically with shooting method

Question

This question is related a question I previously asked here

Solving integro-differential equation with boundary condition at infinity

and for which a solution was found . Now I am dealing with a similar problem with different parameters and I can't get a solution although the asymptotic condition is satisfied. The differential equation contains a complicated integral and has to be solved with a shooting method since we have boundary conditions at infinity:s[0]=0 and s[end]=1. I have 2 different sets of parameters and I'm using the code provided in the topic linked above to get the solution

 (*Set 1*)
\[Chi] = Rationalize[3.5, 0]
A = Rationalize[1.6, 0];
a = Rationalize[1.9123, 0];
b = Rationalize[-28.9815, 0];
d = Rationalize[9.5, 0]; 
n = Rationalize[1, 0]
B = Rationalize[1, 0]


(*Set 2*)
\[Chi] = Rationalize[0.548337, 0]
A = Rationalize[1.1809307701602592, 0];
a = Rationalize[14.593441928706284, 0];
b = Rationalize[-23.834284974146534, 0];
d = Rationalize[5.677752578949083, 0];
n = 0;
(*B is useless*)

(*Set 3*)
χ = Rationalize[0.2, 0];
A = Rationalize[3.915543336008174, 0];
a = Rationalize[181.22232071532673, 0];
b = Rationalize[124.88322492316286, 0];
d = Rationalize[-105.81583286477883, 0];
n = Rationalize[0.06779142517329072, 0];
B = Rationalize[0.479615329820981, 0];
end = 20;

I am using a code slightly modified from the one given by @bbgodfrey in the link given above

 end=12;
 eps=10^-5;
 j[x_, r_] = 2 A^2 r x;
 g0[x_, r_] = a + b/2 + 3 d/4 + A^2 (b + d) (x^2 + r^2) + d A^4 ((x^2 + r^2)^2 + 4 x^2 r^2);
 g1[x_, r_] = j[x, r] (b + 2 d) + 4 d A^4 r x (r^2 + x^2);
 jb[x_, r_] = 2 B^2 r x;
 g0b[x_, r_] = n
 s[0][x_] = x/Sqrt[x^2 + 1/4];
 FNB[0] = Interpolation@Rationalize[Table[{r,2 (Pi)^(3/2)/ A E^(-A^2 (r^2)) NIntegrate[(x s[0][x]^2 E^(-A^2 x^2) (g0[x, r] BesselI[0, j[x, r]] - 
         g1[x, r] BesselI[1, j[x, r]])), {x, 0, 30}] + 2 (Pi)^(3/2)/B E^(-B^2 (r^2)) NIntegrate[(x s[0][x]^2 E^(-B^2 x^2) (g0b[x, r] BesselI[0, jb[x, r]])), {x, 0, 30}]}, {r, 0, end, .1}], 0];

mmin = 1; mmax = 15; imax = 200; wp = 75;

Row[{Dynamic[m], "   ", ProgressIndicator[Dynamic[ip], {0, imax}], "   ", ProgressIndicator[Dynamic[rm], {0, end}]}]

Do[eqnNB = u''[r] + u'[r]/r - u[r]/r^2 + u[r] - \[Chi]*u[r]^5 - 
u[r] (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 == 0;
sp = ParametricNDSolveValue[{eqnNB, u[eps] == 0, u'[eps] == up0, 
WhenEvent[u[r] > 12/10, {bool = 1, "StopIntegration"}], 
WhenEvent[{u[r] < 8/10, u[r] < 0}, {bool = 0,"StopIntegration"}]}, u, {r, eps, end + 1}, {up0, wp0}, 
WorkingPrecision -> wp0, Method -> "StiffnessSwitching",Method -> {"ParametricSensitivity" -> None}, MaxSteps -> 100000];
bl = 1; bu = 10;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
rm = st["Domain"][[1, 2]];
If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i;
If[bool == -1, Return[]], {i, imax}];
s[m] = st; N[bmiddle, wp];
FNB[m] = Interpolation@Rationalize[Table[{r,2 (Pi)^(3/2)/A E^(-A^2 (r^2))      NIntegrate[(x Piecewise[{{s[m][x], 
           eps < x < end}}, 
         s[0][x]]^2 E^(-A^2 (x^2)) (g0[x, r] BesselI[0, j[x, r]] -
          g1[x, r] BesselI[1, j[x, r]])), {x, 0, 30}] + 
   2 (Pi)^(3/2)/
     B E^(-B^2 (r^2)) NIntegrate[(x Piecewise[{{s[m][x], 
           eps < x < end}}, 
         s[0][x]]^2 E^(-B^2 x^2) (g0b[x, r] BesselI[0, 
          jb[x, r]])), {x, 0, 30}]}, {r, 0, end, .1}], 0];, {m,mmin,mmax}]

Plot[Evaluate@Table[s[m][r], {m, mmax - 5, mmax}], {r, eps, end},PlotRange -> All, AxesLabel -> {r, u}, ImageSize -> Large, LabelStyle ->    {Black, Bold, Medium}]
Plot[Evaluate@Table[FNB[m][r], {m, mmax - 5, mmax}], {r, 0, end},PlotRange -> All, AxesLabel -> {r, "FNB"}, ImageSize -> Large,LabelStyle -> {Black, Bold, Medium}]
Plot[Evaluate@Table[s[m][r]^2, {m, mmax - 5, mmax}], {r, 0.01, end}, PlotRange -> All, AxesLabel -> {"r", "n(r)"}, ImageSize -> Large, 
AxesStyle -> Directive[Black, FontSize -> 17,FontFamily -> "TeX Gyre Pagella Math"], LabelStyle -> Directive[Black, FontSize -> 17, FontFamily -> "TeX Gyre Pagella Math"], ImageSize -> {Large}, PlotRange -> All]

For Set 2 integration gets stuck at the first iteration, probably due to divergence, while for Set 1 integration gets stuck at iteration m=4 and i really don't know why: the previous iteration was going on well and fast and were appearing oscillation in the solution that I expect as in the previous case. Here's a picture of how the solution should be (solid line)

Answer 1

Set 1

The code in the question apparently failed to converge, because bl was too large. As explained in the Addendum to my answer to 176883, bl and bu must bracket the final computed value of up0. As it happens, the final value of up0 is approximately 37/100. Setting

bl = 1/5; bu = 4;

meets this requirement with a significant safety margin. (Too much margin can increase runtime by several percent.) With this change, the converged solution after twelve iterations, which requires a few hours, is

Plot[Evaluate@Table[s[i][r], {i, 11, 12}], {r, eps, Min[rm, end]}, PlotRange -> All, 
    AxesLabel -> {r, u}, ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

Incidentally, setting u[eps] == up0 eps in the definition of sp is slightly more accurate than setting u[eps] == 0, because the ODE in the question reduces to the Modified Bessel's Equation of order 1 for very small r. The resulting change to the plotted solution is not visible to the eye but can impact the appropriate value of bl.

Set 2

The code in the question definitely failed for Set 2 because the values of bl and bu did not bracket the final value for up0. A small amount of experimentation indicates that

bl = 3/10; bu = 1;

works well, however. Defining

s[0][x_] = x/Sqrt[x^2 + 1];

instead of the expression in the question accelerates convergence of the iterative solution. Finally, replacing end + 1 by end + 2 in the definition of sp improves computation of the solution near r = end. (Increasing the value of end itself also could be considered.) With these changes, ten iterations gives the converged solution

Plot[Evaluate@Table[s[m][r], {m, 9, 10}], {r, eps, end}, PlotRange -> All, 
    AxesLabel -> {r, u}, ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

as desired.

Addendum

As requested by the OP, the Set 2 results can be computed with end = 20 in a few hours, yielding

Note that NIntegrate has precision issues for larger values of x, which can be circumvented with the options, WorkingPrecision -> 30, PrecisionGoal -> 8.

Set 3

The OP proposed an additional set of parameters in a comment below.

χ = Rationalize[0.2, 0];
A = Rationalize[3.915543336008174, 0];
a = Rationalize[181.22232071532673, 0];
b = Rationalize[124.88322492316286, 0];
d = Rationalize[-105.81583286477883, 0];
n = Rationalize[0.06779142517329072, 0];
B = Rationalize[0.479615329820981, 0];
end = 20;

It appears that the iterative solution was not converging well for r near end. Enough minor changes were needed to resolve this difficulty that reproducing the resulting code here makes sense. With the functions {j, g0, g1, jb, g0b} as defined in the question, the code is

eps = 10^-5; endp = end + 1;
s[0][x_] = x/Sqrt[x^2 + 1];
FNB[0] = Interpolation@Rationalize[Table[{r, 
    2 (Pi)^(3/2)/A E^(-A^2 (r^2)) NIntegrate[(x s[0][x]^2 E^(-A^2 x^2) 
      (g0[x, r] BesselI[0, j[x, r]] - g1[x, r] BesselI[1, j[x, r]])), {x, 0, 30}, 
      WorkingPrecision -> 30, PrecisionGoal -> 8] + 
     2 (Pi)^(3/2)/B E^(-B^2 (r^2)) NIntegrate[(x s[0][x]^2 E^(-B^2 x^2) 
       (g0b[x, r] BesselI[0, jb[x, r]])), {x, 0, 30}, 
      WorkingPrecision -> 30, PrecisionGoal -> 8]}, {r, 0, endp, 1/10}], 0];

mmin = 1; mmax = 10; imax = 200; wp = 60;
Row[{Dynamic[m], "   ", ProgressIndicator[Dynamic[ip], {0, imax}], 
    "   ", ProgressIndicator[Dynamic[rm], {0, end}]}]
Do[eqnNB = u''[r] + u'[r]/r - u[r]/r^2 + u[r] - χ*u[r]^5 - 
        u[r] (FNB[m - 1][r] + FNB[Max[m - 2, 0]][r])/2 == 0;
    sp = ParametricNDSolveValue[{eqnNB, u[eps] == eps up0, u'[eps] == up0, 
        WhenEvent[u[r] > s[0][end] + 1/5 (1 - r/endp) + 1/200 (r/endp), 
            {bool = 1, "StopIntegration"}], 
        WhenEvent[{u[r] < s[0][end] - 1/5 (1 - r/endp) - 1/200 (r/endp), 
            u[r] < 0}, {bool = 0, "StopIntegration"}]}, 
        u, {r, eps, endp}, {up0, wp0}, WorkingPrecision -> wp0, 
        Method -> "StiffnessSwitching", Method -> {"ParametricSensitivity" -> None}, 
        MaxSteps -> 100000];
    bl = 1/5; bu = 2;
    Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp]; 
        rm = st["Domain"][[1, 2]]; 
        If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
        If[bool == -1, Return[]], {i, imax}];
    s[m] = st; N[bmiddle, wp];
    FNB[m] = Interpolation@Rationalize[Table[{r, 
        2 (Pi)^(3/2)/A E^(-A^2 (r^2)) NIntegrate[(x 
            Piecewise[{{s[m][x], eps < x < end}}, s[0][x]]^2 E^(-A^2 (x^2)) 
            (g0[x, r] BesselI[0, j[x, r]] - g1[x, r] BesselI[1, j[x, r]])), 
            {x, 0, 30}, WorkingPrecision -> 30, PrecisionGoal -> 8] + 
        2 (Pi)^(3/2)/B E^(-B^2 (r^2)) NIntegrate[(x 
            Piecewise[{{s[m][x], eps < x < end}}, s[0][x]]^2 E^(-B^2 x^2) 
            (g0b[x, r] BesselI[0, jb[x, r]])), {x, 0, 30}, 
            WorkingPrecision -> 30, PrecisionGoal -> 8]}, {r, 0, endp, 1/10}], 0];,
 {m, mmin, mmax}];

The main changes involve the WhenEvent criteria for stopping integration when a numerical solution overshoots or undershoots the expected asymptotic value of u, the criteria now being much more restrictive at large r. There are a number of lesser changes as well. The well converged solution is

Plot[Evaluate@Table[s[i][r], {i, mmax - 1, mmax}], {r, eps, Min[rm, end]}, 
    PlotRange -> All, AxesLabel -> {r, u}, ImageSize -> Large, 
    LabelStyle -> {Black, Bold, Medium}]