Do-loop with ParametricNDSolveValue not giving expected results

I have written code that using a Do-loop. In the loop I am changing the value of x, v, l and R, and looking at the computed value of tot, which should equal Z. It does not matter what is the values of x, v, l or R are, tot should equal Z. However, the loop gives me a different value of tot. Can anyone help, please?

l1 = 0.81
Z = 500;
x0 = 10;
v0 = 0.02;
ϵ = $MachineEpsilon; l0 = 0.071420.; ps = ParametricNDSolveValue[ {y''[r] + 2 y'[r]/r == -4 π l k Exp[-y[r]], y[ϵ] == y0, y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z l]}, {y, y'}, {r, ϵ, R}, {k}, Method -> "StiffnessSwitching", WorkingPrecision -> 30]; 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]][R], {y0, -10, 0}, Evaluated -> False][[1, 2]]; tot = 4 π nn NIntegrate[r^2 Exp[-First[ps[nn]][r]], {r, 0, R}]; Print[NumberForm[i, 5], " ", NumberForm[tot, 10]];, {i, 2.92, 3.1, 0.01}] • How do we know that Tot = Z? This is not visible from your code. Publish the original model so that we can determine where you are wrong. – Alex Trounev Oct 5 '18 at 4:24 • because this 4 [Pi] nn NIntegrate[r^2 Exp[-First[ps[nn]][r]], {r, 0, R}] is Z – Alz60211 Oct 5 '18 at 4:25 • The code contains errors as reported by the system, for example ParametricNDSolveValue::ndsz: At r == 0.011814321203945158960724620268818668925976400238503304865530., step size is effectively zero; singularity or stiff system suspected. – Alex Trounev Oct 5 '18 at 4:29 • you are right. How to fix it? – Alz60211 Oct 5 '18 at 4:36 • The question is, what do you want to fix? We do not know what problem you solve. – Alex Trounev Oct 5 '18 at 4:43 1 Answer I corrected the code so that in this problem all values should be calculated with a given accuracy. There is an overflow in intermediate calculations, but it apparently does not affect the final result. l1 = 0.81; Z = 500; x0 = 10; v0 = 0.02; \[Epsilon] =$MachineEpsilon;

l0 = 0.0714`20.;

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

2.92 500.05

2.93 500.05

2.94 500.05

2.95 500.05

2.96 500.05

2.97 500.06

2.98 500.06

2.99 500.06

1. 500.06

3.01 500.06

3.02 500.06

3.03 500.06

3.04 500.06

3.05 500.06

3.06 500.06

3.07 500.06

3.08 500.06

3.09 500.05

3.1 500.05

• Thank you so much. – Alz60211 Oct 5 '18 at 14:28