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.0714`20.; 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}]