Why Mathematica is not producing output and taking too much time

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_, \[Theta]_] := r^2 + aa^2 Cos[\[Theta]]^2;  TR[r_, \[Theta]_] := r^2 - 2 M r + aa^2;   gtt[r_, \[Theta]_] := -(1 - (2 M r)/R2[r, \[Theta]]);  gt\[Phi][r_, \[Theta]_] := -(( 2 r M aa Sin[\[Theta]]^2)/    R2[r, \[Theta]]); g\[Phi]\[Phi][   r_, \[Theta]_] := (r^2 +      aa^2 + (2  M r (aa^2) )/       R2[r, \[Theta]] Sin[\[Theta]]^2) Sin[\[Theta]]^2;  grr[r_, \[Theta]_] := R2[r, \[Theta]]/TR[r, \[Theta]];  g\[Theta]\[Theta][r_, \[Theta]_] := R2[r, \[Theta]];   gUtt[r_, \[Theta]_] := -(1/    TR[r, \[Theta]]) (r^2 +      aa^2 + (2  M r (aa^2) )/ R2[r, \[Theta]] Sin[\[Theta]]^2);  gUt\[Phi][r_, \[Theta]_] := -((2 M aa r)/(   TR[r, \[Theta]] R2[r, \[Theta]]));  gU\[Phi]\[Phi][r_, \[Theta]_] := (  TR[r, \[Theta]] - aa^2 Sin[\[Theta]]^2)/(  TR[r, \[Theta]] R2[r, \[Theta]] Sin[\[Theta]]^2);  gUrr[r_, \[Theta]_] := TR[r, \[Theta]]/R2[r, \[Theta]];  gU\[Theta]\[Theta][r_, \[Theta]_] := 1/R2[r, \[Theta]]; M = 1; n = 4; glo = FullSimplify[{ {gtt[r, \[Theta]], 0, 0,       gt\[Phi][r, \[Theta]]}, {0, grr[r, \[Theta]], 0, 0}, {0, 0,       g\[Theta]\[Theta][r, \[Theta]], 0}, {gt\[Phi][r, \[Theta]], 0, 0,       g\[Phi]\[Phi][r, \[Theta]]}}]; gup = FullSimplify[{ {gUtt[r, \[Theta]], 0, 0,       gUt\[Phi][r, \[Theta]]}, {0, gUrr[r, \[Theta]], 0, 0}, {0, 0,       gU\[Theta]\[Theta][r, \[Theta]], 0}, {gUt\[Phi][r, \[Theta]], 0,       0, gU\[Phi]\[Phi][r, \[Theta]]}}];   dglo = Simplify[Det[glo]];  crd = {t, r, \[Theta], \[Phi]};  Xup = {t[\[Tau]], r[\[Tau]], \[Theta][\[Tau]], \[Phi][\[Tau]]}; Vup = {Vt, Vr, V\[Theta], V\[Phi]}; Pup = {Pt[\[Tau]], Pr[\[Tau]], P\[Theta][\[Tau]], P\[Phi][\[Tau]]};  Sup = {{Stt[\[Tau]], Str[\[Tau]], St\[Theta][\[Tau]],      St\[Phi][\[Tau]]},     {Srt[\[Tau]], Srr[\[Tau]], Sr\[Theta][\[Tau]], Sr\[Phi][\[Tau]]},    {S\[Theta]t[\[Tau]], S\[Theta]r[\[Tau]], S\[Theta]\[Theta][\[Tau]],      S\[Theta]\[Phi][\[Tau]]},    {S\[Phi]t[\[Tau]], S\[Phi]r[\[Tau]], S\[Phi]\[Theta][\[Tau]],      S\[Phi]\[Phi][\[Tau]]}};   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]], \[Tau]] == Vup[[a]] , {a, 1, n}];    EOM2 = Table[     D[Pup[[a]], \[Tau]] + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(b = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(c =           1\), \(n\)]christoffel[\([\)\(a, b, c\)\(]\)]*         Pup[\([\)\(b\)\(]\)]*Vup[\([\)\(c\)\(]\)]\)\) == -(1/2) \!\( \*UnderoverscriptBox[\(\[Sum]\), \(b = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(d = 1\), \(n\)]riemann[\([\)\(a,            b, c, d\)\(]\)]*Vup[\([\)\(b\)\(]\)]*          Sup[\([\)\(c, d\)\(]\)]\)\)\),    {a, 1, n}];  EOM3 = Table[     D[Sup[[a, b]], \[Tau]] + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(d =           1\), \(n\)]christoffel[\([\)\(a, c, d\)\(]\)]*         Sup[\([\)\(c, b\)\(]\)]*Vup[\([\)\(d\)\(]\)]\)\) + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(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*\[Mu]^2 + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(l =           1\), \(4\)]\((loriemann[\([\)\(i, j, k,            l\)\(]\)]*\((Sup[\([\)\(i, j\)\(]\)])\)*\ \((Sup[\([\)\(k,             l\)\(]\)])\))\)\)\)\)\);  Wvec = Table[2/(\[Mu]*Wfactor)*(\!\( \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(l = 1\), \(4\)]Sup[\([\)\(j,             i\)\(]\)]*           Pup[\([\)\(k\)\(]\)]*\((loriemann[\([\)\(i, k, l,              m\)\(]\)])\)*\((Sup[\([\)\(l, m\)\(]\)])\)\)\)\)\)), {j,      1, n}];   NN = 1/Sqrt[1 - \!\( \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(k =         1\), \(4\)]\((glo[\([\)\(i, k\)\(]\)])\)*Wvec[\([\)\(i\)\(]\)]*       Wvec[\([\)\(k\)\(]\)]\)\)];   {Vt, Vr, V\[Theta], V\[Phi]} = NN (Wvec + Pup);  EOM = Flatten[    Join[{EOM1, EOM2, EOM3} /.          r -> r[\[Tau]] /. \[Theta] -> \[Theta][\[Tau]] /.        Derivative[1][r[\[Tau]]][\[Tau]] -> Derivative[1][r][\[Tau]] /.       Derivative[1][\[Theta][\[Tau]]][\[Tau]] ->        Derivative[1][\[Theta]][\[Tau]]]];  INT1 = {t[0] == 0,     r[0] == r0, \[Theta][0] == \[Theta]0, \[Phi][0] == 0}; INT2 = {Pt[0] == 1.32288, Pr[0] == 0, P\[Theta][0] == 0,     P\[Phi][0] == 0.07143}; INT3 = {{Stt[0] == 0, Str[0] == 0, St\[Theta][0] == 0,      St\[Phi][0] == 0},     {Srt[0] == 0, Srr[0] == 0, Sr\[Theta][0] == 0, Sr\[Phi][0] == 0},    {S\[Theta]t[0] == 0, S\[Theta]r[0] == 0, S\[Theta]\[Theta][0] == 0,      S\[Theta]\[Phi][0] == 0},    {S\[Phi]t[0] == 0, S\[Phi]r[0] == 0, S\[Phi]\[Theta][0] == 0,      S\[Phi]\[Phi][0] == 0}}; INT = Flatten[Join[{INT1, INT2, INT3}]]; r0 = 7; \[Theta]0 = Pi/2; \[Mu] = 1; aa = 0.5; M = 1;  NDSolve[Flatten[Join[{EOM, INT}]], {t, r, \[Theta], \[Phi], Pt, Pr,    P\[Theta], P\[Phi], Stt, Str, St\[Theta], St\[Phi], Srt, Srr,    Sr\[Theta], Sr\[Phi],   S\[Theta]t, S\[Theta]r, S\[Theta]\[Theta], S\[Theta]\[Phi],    S\[Phi]t, S\[Phi]r, S\[Phi]\[Theta], S\[Phi]\[Phi]}, {\[Tau], 0,    1000}]