# Mathematica gives a solution to a differential equation that does not verify it

The equation I want to solve numerically with WolframCloud is:

z[t_] := Exp[2*Sqrt[L/3]*t]
Eq1 = -40/9*L-20/3*y''[t]/y[t]+2/3*y'[t]^2/y[t]^2-6*z'[t]*y'[t]/(z[t]*y[t])+11/6*z'[t]^2/z[t]^2-360/7*1/y[t]^2+97/21*1/(y[t]^2*z[t])==0


I want to solve it for y[t] only, and then verify that it gives y[t] such that if I want to solve this equation for z[t], it returns the correct exponential function. To do so, I use the following code:

Clear[F,R,y]
L = 0.038;
sol0 = NDSolve[{Eq1, y[0]==0.1, y'[0]==5.5},y[t],{t, 0, 1}];
R[t_]:=Evaluate[y[t]/.sol0]
Eq2 =-40/9*L-20/3*R''[t]/R[t]+2/3*R'[t]^2/R[t]^2-6*F'[t]*R'[t]/(F[t]*R[t])+11/6*F'[t]^2/F[t]^2-360/7*1/R[t]^2+97/21*1/(R[t]^2*F[t])==0;
sol = NDSolve[{Eq2,F[0]==1},F[t],{t,0,10},MaxSteps->10000][[1]];
Plot[{Log[Evaluate[F[t]/.sol]], 2*Sqrt[L/3]*t},{t,0,1}]


As you see, I plot only the log of the function so it should give two superimposed lines. But it gives me this:

We can see that until $$\sim 0.6$$ it gives the correct behavior, but I can't explain what is going on with my code to give a wrong answer... Have I done something wrong in my code?

EDIT: After having played with the parameters, it seems that the point where F[t] deviates from the line $$2 \sqrt{\frac{L}{3}} t$$ is exactly the middle of the domain of definition of R[t]. In fact, R[t] has roughly the shape of $$-at^2+b$$, and so its domain of definition has a maximum right at its half. I don't know what it means but if anyone has an idea...

• Well, one obvious issue is that you solve first ode for up to $t=1$. Then you used the solution in the second ode which you now solved for up to $t=10$. So it will extrapolate the first solution, which is not good. Try to solve first ode for $t=10$. But it is stiff ode. So you can try reducing AccuracyGoal and add Method -> {"StiffnessSwitching"}, to help., Your second ode is also stiff. May be more special options are needed. I get in V 13.0.1 NDSolve::ndsz: At t == 1.3637473398450544, step size is effectively zero; singularity or stiff system suspected. Commented Feb 23, 2022 at 8:04
• @Nasser Thank you for your comment! I've tried to change to $t=10$ in the first equation and to add the Method, but it gives the same plot at the end... Commented Feb 23, 2022 at 8:09
• You can simplify your code a little by writing R = NDSolveValue[{Eq1, y[0] == 0.1, y'[0] == 5.5}, y, {t, 0, 1}] and now you do not need to do R[t_]:=Evaluate[y[t]/.sol0] Commented Feb 23, 2022 at 8:21
• @Nasser I tried to do so, but it gives me errors... Commented Feb 23, 2022 at 8:57
• @Nasser I cannot open the notebook because I didn't subscribe to a plan in WolframCloud, but did you get the right straight lines in the plot? Commented Feb 23, 2022 at 9:19

sol0 = NDSolve[{Eq1, y[0] == 0.1, y'[0] == 5.5}, y[t], {t, 0, 1}, InterpolationOrder -> All];

The difficulty with the original code is that R[t] is a third order interpolation function. Hence, R''[t] is only first order, which apparently is not accurate enough to recover z[t]. The option InterpolationOrder -> All leads to more accurate interpolation of higher derivatives.