# Break NDSolve when condition is “almost” constant

I'm trying the solve a differential equation using NDSolve, starting at x=0.001. I'd like the solver to stop when the function

CSq[x_]:=psi[x]^2*x^(2 l + 2) + psi[x - Pi/2]^2*(x - Pi/2)^(2 l + 2)


is "almost" constant, up to a specified tolerance value. My code so far is

tol = 10^-4;
xmin = 0.0001;
l=1;
V[x_] := -1/x Exp[-x];
eqn = {
psi''[x] + 2 (l + 1) psi'[x]/x + psi[x] - V[x] psi[x] == 0,
psi[xmin] == 1,
psi'[xmin] == 0,
CSq[x_] :=
psi[x]^2*x^(2 l + 2) + psi[x - Pi/2]^2*(x - Pi/2)^(2 l + 2),
WhenEvent[CSq1 - CSq[x] < tol, xmax = x, "StopIntegration"],
CSq1=CSq[x]
};
NDSolve[eqn, psi[x], {x, xmin, \[Infinity]}];


where i want the WhenEvent to trigger when to previous value of psi is close enough to the current value. This doesn't work, and leads me to the following questions:

1. In the error message, it says Equation or list of equations expected instead of Null in the first argument, so i guess that i can't define functions inside NDSolve like this. Is there some solution to this?
2. Even if i remove the function declaration and use psi directly everywhere, i still get error like the one above. I suspect that my last line in eqn is what bothers NDSolve. Is there any convenient way to access the previous value and compare it to the current, during the NDSolve operation?

• What is a typical value for CSq1? – Michael Seifert Feb 24 '17 at 16:55
• @MichaelSeifert Roughly 4.6. – bbgodfrey Feb 24 '17 at 23:18
• When x<Pi/2, how to calculate psi[x]? – xzczd Feb 25 '17 at 8:44
• @xzczd Good question. I worked around this issue by setting CSq to zero there. – bbgodfrey Feb 25 '17 at 8:56
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Feb 25 '17 at 13:16

This question can be solved by means of DiscreteVariables.

tol = 10^-3;
xmin = 10^-4;
l = 1;
V[x_] := -1/x Exp[-x];
eqn = {psi''[x] + 2 (l + 1) psi'[x]/x + psi[x] - V[x] psi[x] == 0,
psi[xmin] == 1, psi'[xmin] == 0,
WhenEvent[Mod[x, Pi/2], {tt = psi[x]^2*x^(2 l + 2),
If[Abs[tem[x] + tt - CSq[x]] < tol, "StopIntegration"],
CSq[x] -> tem[x] + tt, tem[x] -> tt}], CSq[xmin] == 0, tem[xmin] == 0};
s = Flatten@NDSolve[eqn, {psi[x], CSq[x]}, {x, xmin, 100}, DiscreteVariables -> {CSq, tem];


Basically, tem stores old values of psi[x]^2*x^(2 l + 2), and CSq stores old values of the corresonding quantity defined in the question, with new values computed when x is a multiple of Pi/2. The computation stops when the absolute value of the change in CSq is less than tol.

Next, extract the Plot limits. (The final number is the value of x at which the integration stopped.)

plotlim = Join[{x}, Flatten[CSq[x] /. s /. x -> "Domain"]]
(* {x, 0.0001, 61.2611} *)

Plot[CSq[x] /. s, plotlim, PlotRange -> All]


Plot[psi[x] /. s, plotlim, PlotRange -> All]


The corresponding calculation for tol = 10^-4 requires WorkingPrecision -> 30 to obtain adequate accuracy from NDSolve, which then stops at x == 179.070781254618214592370672859.