I am interested in solving a system of Differential Equations, particularly


This system has for initial conditions $u(0)=1/\sqrt{2\cdot k_{55}}$, $u'(0)=-i\cdot k_{55}/\sqrt{2\cdot k_{55}}$. As you may have noticed, $f$ is in terms of $t$, and that's because $\tau$ is a function in terms of $t$.

sol = NDSolve[{x''[t] + 3 x'[t] Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4] + 
Sqrt[3/2] Exp[-x[t] Sqrt[2/3]] (1 - Exp[-x[t] Sqrt[2/3]]) == 0, 
a'[t]/a[t] == Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4],
τ'[t] == 1/a[t],
x'[0] == -0.008226306418212731,
x[0] == 5.630991866033891,
a[0] == 1,
τ[149.4517772937791] == 0},
{x, τ, a},
{t, 0, 500}]
a = a[t] /. sol;
\[Tau] = \[Tau][t] /. sol;
x = x[t] /. sol;
xt = x'[t] /. sol;
att = a''[t] /. sol;
tEnd = t /. FindRoot[(att == 0), {t, 0, 150}]
aEnd = a /. t -> tEnd;
H = Sqrt[(xt)^2/6 + (1 - Exp[-x*Sqrt[2/3]])^2/4];
V = 3/4 (1 - Exp[-Sqrt[2/3] x])^2;
Vt = Sqrt[3/2] Exp[-x Sqrt[2/3]] (1 - Exp[-x Sqrt[2/3]]);
Vtt = -Exp[-2 Sqrt[2/3]*x]*(Exp[Sqrt[2/3]*x] - 2);
t55 = t /. FindRoot[(aEnd/a == E^55), {t, 0, 150}];
\[Tau]55 = \[Tau] /. t -> t55;
a55 = a /. t -> t55;
eps1 = (1/2) (Vt/V)^2;
eta = Vtt/V;
eps2 = -4 eps1 + 2 eta;
k55 = a55*(H /. t -> t55);
\[Nu] = 3/2 + eps1 + 1/2*eps2;
f = k55^2 - (\[Nu]^2 - 1/4)/(\[Tau]^2);
solu = NDSolve[{u''[t] == -u[t]*f, u[t55] == Exp[I*k55*\[Tau]55]/Sqrt[2*k55], u'[t55] 
== -I*k55*Exp[I*k55*\[Tau]55]/(a55*Sqrt[2*k55])}, u[t], {t, t55, tEnd}]

My problem here is that, when trying to get solu (run the NDSolve for $u(t)$), it takes too much time to even give a solution. Next, I show the plot of $f$ from $t=0$ to $t=t_{\text{End}}$

enter image description here

Anyway, I am not interested in a solution on that interval but from $t=t_{55}$ to $t=t_{\text{End}}$, so the initial conditions would be u[t]=Exp[-I*k55*\[Tau]55]/Sqrt[2*k55], and u'(t_{55})=-I*k55*Exp[-I*k55*\[Tau]55]/Sqrt[2*k55]. I hope you can recommend any method to get a solution out of the proposed differential equation, since I need the plot of it. The solution has to look like this:

enter image description here

It only gets to 22 using

solu = NDSolve[{u''[t] == -u[t]*f, u[0] == 1/Sqrt[2*k55], 
u'[0] == -I*k55/(a55*Sqrt[2*k55])}, u[t], {t, 0, 22}, 
Method -> "StiffnessSwitching"]

NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 


solu = NDSolve[{u''[t] == -u[t]*f, u[t55] == Exp[-I*k55*\[Tau]55]/Sqrt[2*k55], u'[t55] == -I*k55*Exp[-I*k55*\[Tau]55]/(a55*Sqrt[2*k55])}, u[t], {t, t55, tEnd}]

Gives the following error: NDSolve::ndsz: At t == 33.42969888667501', step size is effectively zero; singularity or stiff system suspected. And when using MaxSteps -> 10^6 it shows a solution from $t=t_{55}$ to $t = 33.429738218477798903'20$.


I took the anzats $u(t)=e^{\int v(t)dt}$, therefore, I got that $v'+v^2+f=0$. I got $v$ using

solv = NDSolve[{v'[t] + (v[t])^2 + f == 0, v[0] == -I*k55}, v[t], {t, 0, 150}]

However, when doing ur = Integrate[v, t], where v=v[t]/.solv, it does not work and does not give any integral at all, not even ur = NIntegrate[v, {t,0,x}]. How do I solve it to get it? With this I will finally get $u(t)$.

enter image description here

  • 1
    $\begingroup$ It's going to take several steps per oscillation. Looks like lots of oscillations. So maybe MaxSteps -> 10^6 or higher. $\endgroup$
    – Goofy
    Commented Jan 27 at 3:46
  • $\begingroup$ @Goofy Yeah, but the think is that I already tried with MaxSteps -> 10^6, and even with MaxSteps -> 10^9. However, I think it is too much for my laptop, and Mathematica just closes randomly. That's why I would like to know if there is any other method to simplify it. $\endgroup$ Commented Jan 27 at 15:37
  • $\begingroup$ There have been several problems like this (highly oscillatory) posted here, but I don't recall any being solved. You would think such a common problem would have been researched. Maybe try on one of the math sites, scicomp or math or mathoverflow. If there is a math solution, maybe it can be implemented in Mathematica. $\endgroup$
    – Goofy
    Commented Jan 27 at 16:34
  • 1
    $\begingroup$ f as computed in your question, is a complicated expression enclosed in curly brackets, in other words it is a List. As a result, using it in NDSolve produces unpredictable results. Use First[f] to correct this problem. With that correction , you will see that NDSolve[{u''[t] == -u[t]*First[f], u[0] == 1/Sqrt[2*k55], u'[0] == -I*k55/(a55*Sqrt[2*k55])}, u[t], {t, 0, 10^-5}] // Flatten produces a rapidly oscillating solution of approximately uniform amplitude near t = 0, as it should. $\endgroup$
    – bbgodfrey
    Commented Jan 29 at 14:24
  • 2
    $\begingroup$ If you hope to have readers help you, you should first correct known errors in your code. By the way, you question largely duplicates mathematica.stackexchange.com/q/296808/1063. $\endgroup$
    – bbgodfrey
    Commented Jan 29 at 15:23


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.