Why Mathematica is not producing output and taking too much time

Question

I'm trying to solve the given system of ODES but the Mathematica is taking too much time and not producing any output. I was trying to check the error by evaluating one one command but there was no error in any command but the equations EOM2, and EOM3 was taking too much time when I was trying to evaluate the equations.

For simple case aa=0, code works, but when I take non-zero aa, it takes a long time and didn't produce output.

Can anyone please guide me how can I fix this problem? Is there any command in Mathematica that can be used to obtain the fast output?

 R2[r_, θ_] := r^2 + aa^2 Cos[θ]^2; 
TR[r_, θ_] := r^2 - 2 M r + aa^2;

 gtt[r_, θ_] := -(1 - (2 M r)/R2[r, θ]); 
gtϕ[r_, θ_] := -(( 2 r M aa Sin[θ]^2)/ 
  R2[r, θ]);
gϕϕ[
  r_, θ_] := (r^2 + 
    aa^2 + (2  M r (aa^2) )/ 
     R2[r, θ] Sin[θ]^2) Sin[θ]^2; 
grr[r_, θ_] := R2[r, θ]/TR[r, θ]; 
gθθ[r_, θ_] := R2[r, θ];

 gUtt[r_, θ_] := -(1/
   TR[r, θ]) (r^2 + 
    aa^2 + (2  M r (aa^2) )/ R2[r, θ] Sin[θ]^2); 
gUtϕ[r_, θ_] := -((2 M aa r)/(
  TR[r, θ] R2[r, θ]));

gUϕϕ[r_, θ_] := (
 TR[r, θ] - aa^2 Sin[θ]^2)/(
 TR[r, θ] R2[r, θ] Sin[θ]^2); 
gUrr[r_, θ_] := TR[r, θ]/R2[r, θ]; 
gUθθ[r_, θ_] := 1/R2[r, θ];
M = 1; n = 4;
glo = FullSimplify[{ {gtt[r, θ], 0, 0, 
     gtϕ[r, θ]}, {0, grr[r, θ], 0, 0}, {0, 0, 
     gθθ[r, θ], 0}, {gtϕ[r, θ], 0, 0,
      gϕϕ[r, θ]}}];
gup = FullSimplify[{ {gUtt[r, θ], 0, 0, 
     gUtϕ[r, θ]}, {0, gUrr[r, θ], 0, 0}, {0, 0, 
     gUθθ[r, θ], 0}, {gUtϕ[r, θ], 0, 
     0, gUϕϕ[r, θ]}}];


dglo = Simplify[Det[glo]]; 
crd = {t, r, θ, ϕ};

Xup = {t[τ], r[τ], θ[τ], ϕ[τ]};
Vup = {Vt, Vr, Vθ, Vϕ};
Pup = {Pt[τ], Pr[τ], Pθ[τ], Pϕ[τ]};

Sup = {{Stt[τ], Str[τ], Stθ[τ], 
    Stϕ[τ]}, 
   {Srt[τ], Srr[τ], Srθ[τ], Srϕ[τ]},
   {Sθt[τ], Sθr[τ], Sθθ[τ],
     Sθϕ[τ]},
   {Sϕt[τ], Sϕr[τ], Sϕθ[τ], 
    Sϕϕ[τ]}};


christoffel = 
  Table[(1/2)*
    Sum[(gup[[i, s]])*(D[glo[[s, k]], crd[[j]] ] + 
        D[glo[[s, j]], crd[[k]] ] - D[glo[[j, k]], crd[[s]] ]), {s, 1,
       n}], {i, 1, n}, {j, 1, n}, {k, 1, n}] ;


riemann =
  Table[ D[christoffel[[i, j, l]], crd[[k]] ] - 
    D[christoffel[[i, j, k]], crd[[l]] ] + 
    Sum[christoffel[[s, j, l]] christoffel[[i, k, s]] - 
      christoffel[[s, j, k]] christoffel[[i, l, s]],
     {s, 1, n}], {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}] ;


loriemann = 
  Table[Sum[glo[[i, m]]*riemann[[m, j, k, l]], {m, 1, n}], {i, 1, 
    n}, {j, 1, n}, {k, 1, n}, {l, 1, n}] ;

 EOM1 = Table[ D[Xup[[a]], τ] == Vup[[a]] , {a, 1, n}];


 EOM2 = Table[ 
   D[Pup[[a]], τ] + \!\(
\*UnderoverscriptBox[\(∑\), \(b = 1\), \(n\)]\(
\*UnderoverscriptBox[\(∑\), \(c = 
         1\), \(n\)]christoffel[\([\)\(a, b, c\)\(]\)]*
        Pup[\([\)\(b\)\(]\)]*Vup[\([\)\(c\)\(]\)]\)\) == -(1/2) \!\(
\*UnderoverscriptBox[\(∑\), \(b = 1\), \(n\)]\(
\*UnderoverscriptBox[\(∑\), \(c = 1\), \(n\)]\(
\*UnderoverscriptBox[\(∑\), \(d = 1\), \(n\)]riemann[\([\)\(a, 
          b, c, d\)\(]\)]*Vup[\([\)\(b\)\(]\)]*
         Sup[\([\)\(c, d\)\(]\)]\)\)\),
   {a, 1, n}];
 EOM3 = Table[ 
   D[Sup[[a, b]], τ] + \!\(
\*UnderoverscriptBox[\(∑\), \(c = 1\), \(n\)]\(
\*UnderoverscriptBox[\(∑\), \(d = 
         1\), \(n\)]christoffel[\([\)\(a, c, d\)\(]\)]*
        Sup[\([\)\(c, b\)\(]\)]*Vup[\([\)\(d\)\(]\)]\)\) + \!\(
\*UnderoverscriptBox[\(∑\), \(c = 1\), \(n\)]\(
\*UnderoverscriptBox[\(∑\), \(d = 
         1\), \(n\)]christoffel[\([\)\(b, c, d\)\(]\)]*
        Sup[\([\)\(a, c\)\(]\)]*Vup[\([\)\(d\)\(]\)]\)\) == 
    Pup[[a]]*Vup[[b]] - Pup[[b]]*Vup[[a]],
   {a, 1, n}, {b, 1, n}];


 Wfactor = 4*μ^2 + \!\(
\*UnderoverscriptBox[\(∑\), \(i = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(j = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(k = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(l = 
         1\), \(4\)]\((loriemann[\([\)\(i, j, k, 
          l\)\(]\)]*\((Sup[\([\)\(i, j\)\(]\)])\)*\ \((Sup[\([\)\(k, 
           l\)\(]\)])\))\)\)\)\)\);

Wvec = Table[2/(μ*Wfactor)*(\!\(
\*UnderoverscriptBox[\(∑\), \(i = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(k = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(m = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(l = 1\), \(4\)]Sup[\([\)\(j, 
           i\)\(]\)]*
          Pup[\([\)\(k\)\(]\)]*\((loriemann[\([\)\(i, k, l, 
            m\)\(]\)])\)*\((Sup[\([\)\(l, m\)\(]\)])\)\)\)\)\)), {j, 
    1, n}];

 NN = 1/Sqrt[1 - \!\(
\*UnderoverscriptBox[\(∑\), \(i = 1\), \(4\)]\(
\*UnderoverscriptBox[\(∑\), \(k = 
       1\), \(4\)]\((glo[\([\)\(i, k\)\(]\)])\)*Wvec[\([\)\(i\)\(]\)]*
      Wvec[\([\)\(k\)\(]\)]\)\)];


{Vt, Vr, Vθ, Vϕ} = NN (Wvec + Pup);

EOM = Flatten[
   Join[{EOM1, EOM2, EOM3} /. 
        r -> r[τ] /. θ -> θ[τ] /. 
      Derivative[1][r[τ]][τ] -> Derivative[1][r][τ] /. 
     Derivative[1][θ[τ]][τ] -> 
      Derivative[1][θ][τ]]];

INT1 = {t[0] == 0, 
   r[0] == r0, θ[0] == θ0, ϕ[0] == 0};
INT2 = {Pt[0] == 1.32288, Pr[0] == 0, Pθ[0] == 0, 
   Pϕ[0] == 0.07143};
INT3 = {{Stt[0] == 0, Str[0] == 0, Stθ[0] == 0, 
    Stϕ[0] == 0}, 
   {Srt[0] == 0, Srr[0] == 0, Srθ[0] == 0, Srϕ[0] == 0},
   {Sθt[0] == 0, Sθr[0] == 0, Sθθ[0] == 0,
     Sθϕ[0] == 0},
   {Sϕt[0] == 0, Sϕr[0] == 0, Sϕθ[0] == 0, 
    Sϕϕ[0] == 0}};
INT = Flatten[Join[{INT1, INT2, INT3}]];
r0 = 7; θ0 = Pi/2; μ = 1; aa = 0.5; M = 1;

NDSolve[Flatten[Join[{EOM, INT}]], {t, r, θ, ϕ, Pt, Pr, 
  Pθ, Pϕ, Stt, Str, Stθ, Stϕ, Srt, Srr, 
  Srθ, Srϕ,
  Sθt, Sθr, Sθθ, Sθϕ, 
  Sϕt, Sϕr, Sϕθ, Sϕϕ}, {τ, 0, 
  1000}]

Answer 1

As I noted in a comment above, EOM as computed in the question takes over 30 minutes and has a LeafCount of 1177079119. These values can be reduced by FullSimplifying christoffel. riemann. and loriemann. For instance,

christoffel = FullSimplify@Table[...]

I also replaced Simplify by FullSimplify in the definition of dglo, although doing so did not have a significant impact. With these changes EOM can be computed in a few minutes, and

LeafCount[EOM]
(* 23063610 *)

which, will still large, is much smaller than before. Next, for convenience, define

var = Join[Xup, Pup, Flatten[Sup]];

and solve the 24 ODEs.

SetSystemOptions["NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 100.`];
NDSolveValue[Flatten[Join[{EOM, INT}]], var, {τ, 0, 1000}];

SetSystemOptions is needed to keep NDSolve from timing out with the initialization error message, "NDSolve::ntdv", as explained here. A plot of the solutions is given by

Plot[Evaluate[%[[;; 8]]], {τ, 0, 250}, PlotRange -> {Automatic, 10}, ImageSize -> Large,
    PlotLegends -> Placed[ToString /@ var, {.9, .5}], LabelStyle -> {15, Black, Bold}]

The remaining variables are identically zero.

%% /. τ -> 1000
(* {1269.97, 8.32087, 1.5708, 56.0725, 1.2494, -0.0268073, 7.30725*10^-17, 
    0.0500981, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.} *)