I am trying to solve a system of second order differential equation, over time, which has periods of fast changing profile, and I need to evaluate the equation around that transition. I have to run it over a scan of points, thus time per point needs to be below O(1) second. The trouble is that for some of the values, the ODE doesn't wanna solve and spits "Indeterminate" or "Complex Infinity". I tried using MaxStepSize, but that only leads to bogus values as results. I also tried to increase the time interval over which my equation is being evaluated, and although it is giving sensible answers, it takes O(1000) seconds for a single point, making it untenable.
I am attaching the code for one specific point, any suggestions would be greatly appreciated. Let me know if there are any further clarifications needed.
(*VALUES FOR VARIABLES USED*)
fa = 1.13*10^13;
Rf = 0.0104;
\[Theta] = 1;
\[Phi] = 1;
Mpl = 12*10^18;
\[CapitalLambda]qcd = 2/10;
f\[CurlyPhi] = Rf* fa;
ma$t[t_] =
If[t < Mpl/(2 \[CapitalLambda]qcd^2), ma (\[CapitalLambda]qcd/Sqrt[Mpl/(2 t)])^4, ma];
t1 = (Mpl^2/(8 ma \[CapitalLambda]qcd^4))^(1/3);
ma = 57*10^-16 (10^12/fa);
Rm = Rf^4 (3/10 87/10 24/10 1/1000 (fa/10^12)^(3/2))/(11/100);
m\[CurlyPhi] = Rm ma;
(* Mass matrix eigenvalues *)
mH$t[t_] = Sqrt[1/2 ma$t[t]^2 (1 + Rm^2 ma^2/ma$t[t]^2 (1 + Rf^2 +
Sqrt[(1 - Rf^2 - 1/Rm^2 ma$t[t]^2/ma^2)^2 + 4 Rf^2]))];
mL$t[t_] = Sqrt[1/2 ma$t[t]^2 (1 + Rm^2 ma^2/ma$t[t]^2 (1 + Rf^2 -
Sqrt[(1 - Rf^2 - 1/Rm^2 ma$t[t]^2/ma^2)^2 + 4 Rf^2]))];
(*Initial and final solving time, increasing the final times makes the equation take forever*)
ti = t1/1000;
tf = t1 250;
(*DIFFERENTIAL EQUATION*)
eqns = {a''[t] + 3/(2 t) a'[t] + ma$t[t]^2 a[t] + m\[CurlyPhi]^2 Rf^2 a[t] +
m[CurlyPhi]^2 Rf\[CurlyPhi][t] == 0, \[CurlyPhi]''[t] + 3/(2 t) \[CurlyPhi]'[t] +
m\[CurlyPhi]^2 \[CurlyPhi][t] + m\[CurlyPhi]^2 Rf a[t] == 0, a'[ti] == 0,
\[CurlyPhi]'[ti] == 0, a[ti] == \[Theta] fa, \[CurlyPhi][ti] ==
\[Phi]f[CurlyPhi]};
{{sol1, sol2}} = NDSolve[eqns, {a[t], \[CurlyPhi][t]}, {t, ti, tf}, MaxSteps -> Infinity,
PrecisionGoal -> 6];
UPDATE : Thanks a lot to Alex, and following his solution involving dimensionless quantities and NDSolve method tricks, it speeded up the evaluation, and also realizing that I had to solve until t0 at max, I used "ADAMS" method. But for some points, I realized that using "BDF" actually resolved the issue better, so would like to understand why's that.
The query can be thought in a generic way : which Methods and tactics are best used for solving 2nd order ODEs, which have solutions which start at some initial value, slowly roll towards it's minima and then oscillate around it.
As an example, am putting in the updated code:
\[Theta] = 1;
\[Phi] = 1;
fa = 2.1544346900318835855901063672841`10.*^12;
ma = 57 10^-16 (10^12/fa); (* GeV axion mass *)
t1 = (Mpl^2/(
8 ma \[CapitalLambda]qcd^4))^(
1/3); (* GeV^-1 time when ma(T) = H(T) *)
t0 = Mpl/(2 \[CapitalLambda]qcd^2)/t1;
Rf = 0.6256495158757365`
fphi = f\[CurlyPhi] = Rf* fa;
Rm =0.0116328;
mphi = m\[CurlyPhi] = Rm ma;
m0[t_] := ma (\[CapitalLambda]qcd/Sqrt[Mpl/(2 t1 t)])^4;
madt[t_] := ma*If[t < t0, m0[t]/ma, 1];
mHDt[t_] :=
Sqrt[1/2 madt[t]^2 (1 +
Rm^2 ma^2/
madt[t]^2 (1 + Rf^2 +
Sqrt[(1 - Rf^2 - 1/Rm^2 madt[t]^2/ma^2)^2 + 4 Rf^2]))];
mLDt[t_] :=
Sqrt[1/2 madt[t]^2 (1 +
Rm^2 ma^2/
madt[t]^2 (1 + Rf^2 -
Sqrt[(1 - Rf^2 - 1/Rm^2 madt[t]^2/ma^2)^2 + 4 Rf^2]))];
ti = 1/1000;
tlc = t0 (1 - Rf^2)^(1/4) Sqrt[
Rm];
tf = 4 tlc; a0 = \[Theta]*fa; phi0 = \[Phi]*fphi;
time = {t, ti, tf};
qvars = {a[t], \[CurlyPhi][t]};
(*DIFFERENTIAL EQUATION*)
eqns = {(a''[t] + 3/(2 t) a'[t]) + t1^2 madt[t]^2 a[t] +
t1^2 mphi^2 Rf^2 a[t] + t1^2 mphi^2 Rf phi0/a0 \[CurlyPhi][t] ==
0, (\[CurlyPhi]''[t] + 3/(2 t) \[CurlyPhi]'[t]) +
t1^2 mphi^2 \[CurlyPhi][t] + t1^2 mphi^2 Rf a0/phi0 a[t] == 0,
a'[ti] == 0, \[CurlyPhi]'[ti] == 0,
a[ti] == 1, \[CurlyPhi][ti] == 1} // Simplify;
solA = NDSolve[eqns, {a, \[CurlyPhi]}, time, StartingStepSize -> 1/20,
Method -> "BDF", InterpolationOrder -> 2,
MaxSteps -> Infinity]; // AbsoluteTiming
```