# Need help solving ODE with ParametricNDSolveValue solver

I am trying to apply @bbgodfrey method for finding the right initial guess for the following problem but I get an error message and I don't understand where it comes from.

I would also like to know why FindRoot in Do loop doesn't work in this case.

Below is my minimal working code:

l1 = 0.81;
Z = 500;
x0 = 10;
v0 = 0.02;
\[Epsilon] = $MachineEpsilon; yl = -12; yu = 0; l0 = 0.071420.; ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == -4 \[Pi] l k Exp[-y[r]], y[\[Epsilon]] == y0, y'[\[Epsilon]] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z l]}, {y, y'}, {r, \[Epsilon], R}, {k, l}, Method -> {"StiffnessSwitching"}, AccuracyGoal -> 5, PrecisionGoal -> 4, WorkingPrecision -> 15]; Do[x = i x0; v = i^3 v0; yl = -12; yu = 0; R = Rationalize[v^(-1/3), 0]; l = Rationalize[l1/(i x0), 0]; fy := (Do[yc = (yl + yu)/2; (* guess finder function *) test = First[ps[yc]]["Domain"][[1, 2]]; If[test == 1, Throw[yc]]; If[Last[ps[yc]][test] > 0, yu = yc, yl = yc], {i, 50}]; yc); yint = Which[1 == First[ps[yl]]["Domain"][[1, 2]], yl, 1 == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]]; nn = FindRoot[Last[ps[y0]][R], {y0, yint}, Evaluated -> False][[1, 2]]; Tot = 4 \[Pi] nn NIntegrate[ r^2 Exp[-First[ps[nn, l]][r]], {r, \[Epsilon], R}, PrecisionGoal -> 4]; Print[NumberForm[i*1., 5], " ", NumberForm[Tot, 5]];, {i, 292/100, 31/10, 1/100}]  • The first error message your code gives me says that ps was called with one parameter having a value of -12. That would be in the line ψint = ... ps[ yl ] ... Your definition of ps expects 2 parameters. However, the definition of ps contains y0 and R, which are undefined. Perhaps y0 and R should be a third and fourth parameters of ps ? If they are not parameters, they need numeric values. Also, to debug loops, start with With[ {i=3}, ... ] and add statements one at a time, making sure each one works for the value of i. Then try to replace the With with Table instead of Do. – LouisB Jan 25 at 3:25 • @LouisB First, I need to work with Do. And in Do statements R is numerical value. For y0, I belived that the code should take yint and pass it to y0. Also, there is a typo ψint should be yint. – aluuzz Jan 25 at 3:37 • Just to understand the syntax, try this immediately after your ps =  statement, y0 = 1; R = 1.2; Plot[ ps[0.5, 0.027] [[1]] [r], {r, 0, R}, PlotRange -> {{0, R}, {0, 4}}]. This indicates we must set values of y0 and R, then evaluate ps with 2 parameters, followed by a subscript of 1 or 2 followed by [ r ]. Of course we could use First and Last instead of the [[ ... ]] notation. – LouisB Jan 25 at 5:00 • @LouisB OK. I just want to let know that the code that I am using is taken from here [mathematica.stackexchange.com/questions/183148/…, and it works fine. What I am asking here is to implemented bbgodfrey method for finding the right initial guess to the above code. – aluuzz Jan 25 at 5:12 • @aluuzz 1) what problem do you want to solve? 2) it looks so that you have syntax errors there. – Alex Trounev Jan 28 at 20:53 ## 1 Answer We can use my code without changing for Z<=15000 l1 = 0.81; Z = 15000; x0 = 10; v0 = 0.02; \[Epsilon] =$MachineEpsilon;

l0 = 0.071420.;

ps = ParametricNDSolveValue[{y''[r] +
2 y'[r]/r == -4 \[Pi] l k Exp[-y[r]], y[\[Epsilon]] == y0,
y'[\[Epsilon]] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z l]}, {y,
y'}, {r, \[Epsilon], R}, {k, l},
Method -> {"StiffnessSwitching"}, AccuracyGoal -> 5,
PrecisionGoal -> 4, WorkingPrecision -> 15];

Do[x = i x0;
v = i^3 v0;
R = Rationalize[v^(-1/3), 0];
l = Rationalize[l1/(i x0), 0];
nn = FindRoot[Last[ps[y0, l]][R], {y0, -1}, Evaluated -> False][[1,
2]];
Tot = 4 \[Pi] nn NIntegrate[
r^2 Exp[-First[ps[nn, l]][r]], {r, \[Epsilon], R},
PrecisionGoal -> 4];
Print[NumberForm[i*1., 5], "  ", NumberForm[Tot, 5]];, {i, 292/100,
31/10, 1/100}] // Quiet


Out:

2.92  15007.

2.93  15007.

2.94  15006.

2.95  15006.

2.96  15006.

2.97  15006.

2.98  15006.

2.99  15006.

3.  15006.

3.01  15006.

3.02  15006.

3.03  15006.

3.04  15006.

3.05  15006.

3.06  15006.

3.07  15006.

3.08  15005.

3.09  15005.

3.1  15005.


I checked code @bbgodfrey (question as it is). The result was much worse than with my code, but the code works for some k up to Z = 20000:

p[Z0_, g0_, k0_, R0_] :=
Block[{Z = Z0, g = Rationalize[g0, 0], k2 = Rationalize[k0, 0],
yl = -8, yu = 0, ps, fy, y00, sol}, \[Epsilon] =
10^-4; R = Rationalize[R0, 0];
ps = ParametricNDSolveValue[{y''[r] +
2 y'[r]/r == -4 Pi k2 Exp[-y[r]], y[\[Epsilon]] == y0,
y'[\[Epsilon]] == 0,
WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y,
y'}, {r, \[Epsilon], R}, {y0}, Method -> "StiffnessSwitching",
WorkingPrecision -> 20];
fy := (Do[yc = (yl + yu)/2;
tst = First[ps[yc]]["Domain"][[1, 2]];
If[tst == R, Throw[yc]];
If[Last[ps[yc]][tst] > 0, yu = yc, yl = yc], {i, 50}]; yc);
y00 = Which[R == First[ps[yl]]["Domain"][[1, 2]], yl,
R == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]];
sol = FindRoot[Last[ps[y0]][R], {y0, y00}, Evaluated -> False][[1,
2]]; L =
4 \[Pi] k2/g0 NIntegrate[
r^2 Exp[-First[ps[sol]][r]], {r, \[Epsilon], R}]]


Using of p[]

l1 = 0.81;
Z0 = 20000;
x0 = 10;
v0 = 0.02;
l0 = 0.071420.; k = 1/29;
Do[x = i x0;
v = i^3 v0;
R = Rationalize[v^(-1/3), 0];
l = Rationalize[l1/(i x0), 0];
Print[NumberForm[i*1., 5], "  ",
NumberForm[p[Z0, l, k l, R] // Quiet, 5]];, {i, 292/100, 31/10,
1/100}]


Out

2.92  20010.

2.93  20010.

2.94  20010.

2.95  20009.

2.96  20009.

2.97  20009.

2.98  20009.

2.99  20008.

3.  20008.

3.01  8.3397*10^(21)

3.02  1.309*10^(22)

3.03  2.0575*10^(22)

3.04  3.2386*10^(22)

3.05  5.1048*10^(22)

3.06  8.0576*10^(22)

3.07  1.2736*10^(23)

3.08  2.0159*10^(23)

3.09  3.1952*10^(23)

3.1  5.0714*10^(23)
`