How to run NDSolve over a large range without freezing?

Question

I am trying to solve a set of coupled differential eqns. I need the solution over a large range of the variable t. But the notebook freezes and stops working.

replace = {Subscript[m, ϕ] -> 10^-5, Γ -> 
    10^-11, λ -> 0.01, ξ -> -1, m -> 10^-21, 
   Subscript[M, P] -> 1};

V[t_] := 1/2 Subscript[m, ϕ]^2  ϕ[t]^2  /. replace ;

R[t_] :=   Subscript[M, P]^-2 (4 V[t] - ϕ'[t]^2) /. replace;

H[t_] := Sqrt[
   1/(3 Subscript[M, 
     P]^2) (1/2 ϕ'[t] ^2  + V[t] + Subscript[ρ, r][t])] /. 
   replace;

χi[t_] := ((- ξ R[t] - m^2)/λ)^(1/2) /. replace;

eqna = ϕ''[t]  + 
    3 H[t] ϕ'[t] + Γ ϕ'[t] + 
    D[V[t], ϕ[t]] /. replace;

eqnb = χ''[t]  + 3 H[t] χ'[t] + λ χ[t]^3 + 
    m^2 χ[t] + ξ R[t] χ[t] /. replace;

eqnc = Subscript[ρ, r]'[t] + 
    4 H[t] Subscript[ρ, r][t] - Γ ϕ'[t]^2 /. 
   replace;
sol1 = NDSolve[{eqna == 0 , eqnb == 0, 
    eqnc == 0,  ϕ[10^5] == 15 , ϕ'[10^5] == 
     0, χ[10^5] == χi[10^5], χ'[10^5] == 0, 
    Subscript[ρ, r][10^5] == 10^-20}, {ϕ, χ, 
    Subscript[ρ, r]}, {t, 10^5, 10^50}] // FullSimplify

The code works fine for smaller intervals: say around 10^5 to 10^10. Is there a way I could run the program in a short time for my required large interval without freezing the notebook? Taking some time, say around 1 day is also fine, provided the notebook doesn't freeze.

Note that I have also tried redefining the variable in Log scale, so that the interval is reduced. But the problem still remains.

replace = {Subscript[m, ϕ] -> 10^-5, Γ -> 
    10^-11, λ -> 0.01, ξ -> -1, m -> 10^-21, 
   k -> Log[10], Subscript[M, P] -> 1};

V[t_] := 1/2 Subscript[m, ϕ]^2  ϕ[t]^2  /. replace ;

R[t_] :=   
  Subscript[M, P]^-2 (4 V[t] - k^-2 10^(-2 t) ϕ'[t]^2) /. replace;

H[t_] := Sqrt[
   1/(3 Subscript[M, 
     P]^2) (1/2 ϕ'[t] ^2  + k^2 10^(2 t) V[t] + 
      k^2 10^(2 t) Subscript[ρ, r][t])] /. replace;

χi[t_] := ((- ξ R[t] - m^2)/λ)^(1/2) /. replace;

eqna = k^-2 10^(-2 t) ϕ''[t] - ϕ'[t] k^-1 10^(-2 t)   + 
    3 k^-2 10^(-2 t)
      H[t] ϕ'[t] + Γ k^-1 10^-t ϕ'[t] + 
    D[V[t], ϕ[t]] /. replace;

eqnb =  k^-2 10^(-2 t) χ''[t]  - χ'[t] k^-1 10^(-2 t)  + 
    3 k^-2 10^(-2 t) H[t] χ'[t] + λ χ[t]^3 + 
    m^2 χ[t] + ξ R[t] χ[t] /. replace;

eqnc = k^-1 10^-t Subscript[ρ, r]'[t] + 
    4 k^-1 10^-t H[t] Subscript[ρ, r][
      t] - Γ k^-2 10^(-2 t) ϕ'[t]^2 /. replace;
sol1 = NDSolve[{eqna == 0 , eqnb == 0, 
    eqnc == 0,  ϕ[5] == 15 , ϕ'[5] == 
     0, χ[5] == χi[5], χ'[5] == 0, 
    Subscript[ρ, r][5] == 10^-20}, {ϕ, χ, 
    Subscript[ρ, r]}, {t, 5, 50}] // FullSimplify 
```

Answer 1

Update: Complete solution for ϕ and ρ

Although my original answer, that Mathematica appears to freeze, because it is running out of memory, is correct, I was troubled by the abrupt changes in the magniturde of ϕ in the first plot of the answer. It turns out that Mathematica's default working precision is inadequate for t > 10^9. The usual approach in such cases, using a larger working precision (say, 30) is not acceptable in this case, because doing so consumes much more memory, the quantity already in short supply. Instead, it is possible to rescale t to mf t (mf replacing Subscript[m, ϕ] for convenience), ρ (replacing Subscript[ρ, r]) by ρ/mf^2, and H by H/mf. Moreover, χ and eqnb can be dropped from the system for now, because χ does not enter the other equations. The system in the question then reduces to

replace = {mf -> 10^-5, Γ -> 10^-11, mp -> 1};
H[t_] := Sqrt[1/3 (1/2 ϕ'[t]^2 + 1/2  ϕ[t]^2 + ρ[t])]/mp
eqna = ϕ''[t] + 3 H[t] ϕ'[t] + Γ/mf ϕ'[t] + ϕ[t]
eqnc = ρ'[t] + 4 H[t] ρ[t] - Γ/mf ϕ'[t]^2
sol3 = NDSolveValue[{eqna == 0, eqnc == 0, ϕ[1] == 15, ϕ'[1] == 0, 
    ρ[1] == 10^-20/mf^2} /. replace, {ϕ, ρ, ϕ'}, {t, 1, 10^7}, 
    AccuracyGoal -> $MachinePrecision, PrecisionGoal -> 13, MaxSteps -> Infinity];

(where mp replaces Subscript[M, P]. In general, it is better not to use subscripted variables in computations.) The transformation and better choice of NDSolve options, allows the computation to proceed accurately to t = 10^7 (10^12 before the transformation`).

Plot[sol3[[1]][t], {t, 1, 10^7}, PlotRange -> {-10^-7, 10^-7}, 
    AxesLabel -> {t, ϕ}, PlotPoints -> 10^3, ImageSize -> Large, 
    LabelStyle -> {Bold, Black, 15}]

LogLogPlot[sol3[[2]][t], {t, 1, 10^7}, PlotRange -> All, 
AxesLabel -> {"t", ρ}, ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

but little further. For very large t, however, it is possible to average over the rapid oscillations in ϕ, which has period 2 Pi in the transformed time coordinate. Doing so reduces the system to

Has[t_] := Sqrt[1/3 (1/2  ϕ[t]^2 + ρ[t])]/mp
eqnas = 2 ϕ'[t] + 3 Has[t] ϕ[t] + Γ/mf ϕ[t]
eqncs = ρ'[t] + 4 Has[t] ρ[t] - Γ/mf ϕ[t]^2/2
sol4 = NDSolveValue[{eqnas == 0, eqncs == 0, ϕ[10^4] == 0.000162926, 
    ρ[10^4] == 7.98079 10^-11} /. replace, {ϕ, ρ}, {t, 10^4, 10^7}];

with starting conditions taken from the solution to sol3. Comparision with the solution of sol3 is very good. For much larger integration domains, lϕ = Log[ϕ] is a better behaved dependent variable.

Hal[t_] := Sqrt[1/3 (1/2 Exp[2 lϕ[t]] + ρ[t])]/mp
eqnal = 2 lϕ'[t] + 3 Hal[t] + Γ/mf
eqncl = ρ'[t] + 4 Hal[t] ρ[t] - Γ/mf Exp[2 lϕ[t]]/2
sol6 = NDSolveValue[{eqnal == 0, eqncl == 0, 
    lϕ[10^4] == Log[0.000162926], ρ[10^4] == 7.98079 10^-11} /. replace, 
    {lϕ, ρ}, {t, 10^4, 10^10}];)

LogLogPlot[Exp[sol6[[1]][t]], {t, 10^4, 10^10}, 
    PlotRange -> {Full, {10^-100, 10^-3}}, AxesLabel -> {"t", ϕ}, 
    ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

LogLogPlot[sol6[[2]][t], {t, 10^4, 10^10}, PlotRange -> All, 
    AxesLabel -> {"t", ρ}, ImageSize -> Large, 
    LabelStyle -> {Bold, Black, 15}]

Finally, for such tiny values of ϕ, the two ODEs decouple and can be solved symbolically.

DSolveValue[Γ/mf ϕ[t] + 2 ϕ'[t] == 0, ϕ[t], t]
(* E^(-((t Γ)/(2 mf))) C[1] *)
DSolveValue[(4 ρ[t]^(3/2))/(Sqrt[3] mp) + ρ'[t] == 0, ρ[t], t]
(* (12 mp^2)/(16 t^2 - 8 Sqrt[3] mp t C[2] + 3 mp^2 C[2]^2) *)

which connect smoothly with the curves of the third and fourth plots with appropriate choice of constants. χ can now be obtained, if desired.

Original Answer

In cases like this, it often helps to plot a partial solution, in this case for {t, 10^5, 10^10}. (I also use NDSolveValue instead of NDSolve for notational simplicity and omit FullSimplify, which has no effect here.)

Plot[sol1[[1]][t], {t, 10^5, 10^10}, PlotRange -> {-10^-5, 10^-5}, AxesLabel -> {t, ϕ}]

which displays oscillations with period of approximately 2π 10^5, which are not properly resolved in this plot.

Plot[sol1[[2]][t], {t, 10^5, 10^10}, PlotRange -> {-10^-6, 10^-6}, AxesLabel -> {t, χ}]

which displays oscillations with period of approximately 10^8.

LogLogPlot[sol1[[3]][t], {t, 10^5, 10^10}, PlotRange -> All, AxesLabel -> {t, Subscript[ρ, r]}]

It is clear from the first two plots that NDSolve requires very many steps in t to represent the solution. To be precise,

sol1[[1]]["Grid"] // Length
(* 229910 *)

It is reasonable to suppose that an upper bound of 10^50 instead of 10^10 would required 10^40 times as many steps, far beyond the memory of any PC. When Mathematica runs out of memory, it often moves the solution to disk, which is excruciatingly slow and appears to lock up the PC.

Computing to 10^10 may be sufficient, because the solution appears to settle into a pattern by then. It also may be possible to find a symbolic solution in the asymptotic regime.