# 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 – aluuzz 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? – aluuzz 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.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
`

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. – aluuzz Oct 5 '18 at 14:28