# How to keep the numerical integration when an event occurs?

I'm numerically integrating the following system,

f[r_] = 0.904965 - (1 + 12.7028/r^2) (1 - 2/r)
dϕdλ = 1/r[λ]^2;
eqc = r'[λ] == Sqrt[f[r[λ]]/(dϕdλ^2)];


with the numerical routine,

sc = NDSolve[{eqc, r[0] == 4.85744 + .003},r, {λ, 0., 40*Pi}, AccuracyGoal -> 10, PrecisionGoal -> 20, Method -> {"StiffnessSwitching"}]


However, I need to keep the integration when the solution reaches the radius r = 4.85744 or r = 11.3301, which are the values where the derivative changes the signal.

Technically, what I need is: When r[λ] = 4.85744 or r[λ] = 11.3301 keep the integration with r'[λ]==-Sqrt[f[r[λ]]/(dϕdλ^2)].

In my code, the integration just stops at this point. I would like to keep integrating with opposite velocity when the solution reaches these two radii, is it possible?

• There was a little mistake but I've fixed and edited the post, I hope I made it clear @ChrisK. – Herr Schrödinger May 14 at 20:14
• Thanks for the edit! – Chris K May 14 at 20:15
• There is a solution with reflection from point r=11.3301, but there is no solution with reflection from point r=4.85744. – Alex Trounev May 14 at 23:25
• @AlexTrounev That’s exactly the solution that I want to integrate! – Herr Schrödinger May 14 at 23:43

It is necessary to square the equation eqc and differentiate it once by lambda. Then we can use WhenEvent[r[\[Lambda]] == 4.85744, r'[\[Lambda]] -> -r'[\[Lambda]]].

F[r_] := 2 (0.904965 - (1 + 12.7028/r^2) (1 - 2/r)) r + ((
25.4056 (1 - 2/r))/r^3 - (2 (1 + 12.7028/r^2))/r^2) r^2
mu[r_] := Sqrt[0.904965 - (1 + 12.7028/r^2) (1 - 2/r)]*r;
eqc = 2*r''[\[Lambda]] == F[r[\[Lambda]]]; r0 = 4.86;
sc = NDSolve[{eqc, r[0] == r0, r'[0] == mu[r0],
WhenEvent[r[\[Lambda]] == 4.85744,
r'[\[Lambda]] -> -r'[\[Lambda]]]}, r, {\[Lambda], 0., 80*Pi},
Method -> {"StiffnessSwitching"}]

Plot[Evaluate[r[\[Lambda]] /. sc], {\[Lambda], 0., 80*Pi},
PlotRange -> {0, 12}, AxesLabel -> {"\[Lambda]", "r"}]