How to directly input calculated Geodesic differential equation into *NDSolve*

Question

I am trying to find the trajectory of a test particle around a relativistic object. Initially, I was trying the Kerr metric as spacetime but the final simulations clearly were not correct as I saw the particle being repelled by the star/planet/black-hole which is not realistic. Hence there must be something wrong going on with the equations I have calculated and used. Since the equations were humongous I was not able to rectify the mistake (although it was not me but my computer that calculated the differential equations, so the equations would be wrongs seems not so plausible to me). I did modify a little bit of the equation I was calculated. I'll show how. (I'll use the Schwarzschild metric and not Kerr metric as the equations are very ugly)

First I defined the metric tensor for the schwarzchild geometry.

schwar = {{-(1 - (2 G M)/(r c^2))^-1, 0, 0, 0}, {0, -r^2, 0, 0}, {0, 
          0, -r^2 Sin[\[Theta]]^2, 0}, {0, 0, 0, (1 - (2 G M )/(c^2 r))}};

then I defined the inverse metric

inverseschwar = Inverse[schwar];

Now, I have to calculate the Christoffel symbols for the metric. I defined a function christ for this

christ[ a_, b_, c_] := christ[a, b, c] =   Sum[1/2 inverseschwar[[a, d]] 
                        (D[schwar[[d, c]], coord[[b]]] + D[schwar[[d, b]], 
                    coord[[c]]] - D[schwar[[b, c]], coord[[d]]]),
                    {d, 1, 4}]

Now I have to compute the geodesic differential equations from the geodesic equation

Manipulate[coord[[i]]'' + Sum[christs[i, h, k] coord[[h]]'
           coord[[k]]', {h, 1, 4}, {k, 1, 4}], {i, {1, 2, 3, 4}}]

This above are the 4 differential equations for each coordinate $r$,$\theta$,$\phi$ and $t$. Notice that I have purposely put on the differentiation operators $''$ and $'$ . But these operators don't do differentiation in mathematica (Once you run this code you'll understand why. They come on on the variable making in into a new variable).

For example: output for the equation of $\theta$ coordinate which is

Here the $\theta$ and $\theta^{'}$ are being considered as two different coordinates and not as differentiation (w.r.t proper time $s$) of $\theta$. Moreover what I need is that all coordinates to be labeled as a function of $s$ so that integration process can be carried out. e.g. $\theta$ should be represented as $\theta[s]$. I cannot do this by just putting an $s$ in front of every variable before calculating these equations. I can only do it for derivative terms present in the equations but cannot for the terms that arise from the chirstoffel symbols as if suppose I put the argument $s$ before calculating the chirstoffel symbols then as the metric is being differentiated w.r.t coordinates $r$, $\theta$, $\phi$ and $t$ so every chirstoffel symbol will become zero as $s$ is independent of all coordinates.

So now I have to manually modify the equationsto make it look like this

After making it look like above now I can feed this and other $3$ equation to NDSolve function and get the output trajectory of the particle.

So, my question is, is there a way to directly have mathematica produce the equation that can directly be inserted into NDSolve function? I don't want to make the process manual at any step.

I am listing the 4 differential equations for your reference

Radial Equation

(G M r'[s]^2)/(2 G M r[s] - c^2 r[s]^2) + (
G M (-2 G M + c^2 r[s]) t'[s]^2)/(
c^4 r[s]^3) + ((2 G M)/c^2 - r[s]) \[Theta]'[
 s]^2 + ((2 G M - c^2 r[s]) Sin[\[Theta][s]]^2 \[Phi]'[s]^2)/c^2 +
r''[s]==0;

Theta Equation

(2 r'[s] \[Theta]'[s])/r[s] - 
Cos[\[Theta][s]] Sin[\[Theta][s]] \[Phi]'[s]^2 + \[Theta]''[s]==0;

Phi Equation

(2 r'[s] \[Phi]'[s])/r[s] + 
2 Cot[\[Theta][s]] \[Theta]'[s] \[Phi]'[s] + \[Phi]''[s]==0;

Time Equation

(2 G M r'[s] t'[s])/(r[s] (-2 G M + c^2 r[s])) + t''[s]==0;

Answer 1

I think the following works (check that it is correct please).

schwar = {{-(1 - (2 G M)/(r c^2))^-1, 0, 0, 0}, {0, -r^2, 0, 0}, {0, 
    0, -r^2 Sin[\[Theta]]^2, 0}, {0, 0, 0, (1 - (2 G M)/(c^2 r))}};   
inverseschwar = Inverse[schwar];
christ[a_, b_, c_] := 
 christ[a, b, c] = 
  Sum[1/2 inverseschwar[[a, d]] (D[schwar[[d, c]], coord[[b]]] + 
      D[schwar[[d, b]], coord[[c]]] - 
      D[schwar[[b, c]], coord[[d]]]), {d, 1, 4}]

Then

christ[1, 1, 1]

coord = {r, \[Theta], \[Phi], t}

Table[D[coord[[i]][s], s, s] + 
     Sum[christ[i, h, k] D[coord[[h]][s], s] D[coord[[k]][s], s], {h, 
       1, 4}, {k, 1, 4}] == 0, {i, {1, 2, 3, 4}}] // 
  FullSimplify // TableForm

should provide your result:

Update

If you want the functions to appear explicitly in the solution

coord = {r, \[Theta], \[Phi], t};
rc = Thread[coord -> Map[#[s] &, coord]];
Clear[christ]
christ[a_, b_, c_] := 
 Sum[1/2 inverseschwar[[a, d]] (D[schwar[[d, c]], coord[[b]]] + 
      D[schwar[[d, b]], coord[[c]]] - 
      D[schwar[[b, c]], coord[[d]]]), {d, 1, 4}] /. rc
christ[1, 1, 1]

Table[D[coord[[i]][s], s, s] + 
    Sum[christ[i, h, k] D[coord[[h]][s], s] D[coord[[k]][s], s], {h, 
      1, 4}, {k, 1, 4}] == 0, {i, {1, 2, 3, 4}}] // FullSimplify