ClearAll["Global`*"] (* geodesics.yukterez.net *)

(* Input: kovariante metrische Komponenten *)
gtt=1-(2r-(℧^2))/Σ;
grr=-Σ/Δ;
gθθ=-Σ;
gφφ=-Χ/Σ ((Sin[θ])^2);
gtφ=a (2r-(℧^2)) ((Sin[θ])^2)/Σ;

(* Abkürzungen *)
Σ=(r^2)+(a^2) ((Cos[θ])^2);
Δ=(r^2)-2 r+(a^2)+(℧^2);
Χ=((((r^2)+(a^2)))^2)-(a^2) ((Sin[θ])^2) Δ;
σ=((r[τ])^2)+(a^2) ((Cos[θ[τ]])^2);
δ=((r[τ])^2)-2 r[τ]+(a^2)+(℧^2);
щ=(q ℧ r[τ] ((a^2)+((r[τ])^2)))/(δ σ);

(* Koordinaten, Dimensionen, magnetisches Monopol, elektrische Ladung, Spin *)
x={t, r, θ, φ}; n=4; P=0; Ω=℧; ℧=Q; a=a;

"Metrischer Tensor"
metrik={{gtt, 0, 0, gtφ}, {0, grr, 0, 0}, {0, 0, gθθ, 0}, {gtφ, 0, 0, gφφ}};
g1=Subscript["g", μσ] -> MatrixForm[metrik];
inversemetrik=Simplify[Inverse[metrik]];g2="g"^μσ -> MatrixForm[inversemetrik];
Grid[{{g1,"  ", g2}}]

"Maxwell Tensor"
A={(Ω r)/Σ+P/Σ Cos[θ] a, 0, 0, -(Ω r/Σ) ((Sin[θ])^2) a-P/Σ Cos[θ]((r^2)+(a^2))};
F=Simplify[Table[((D[A[[j]], x[[k]]]-D[A[[k]], x[[j]]])), {j, 1, 4}, {k, 1, 4}], Reals]; F1=Subscript["F", μσ] -> MatrixForm[F];
f=FullSimplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] F[[k, l]], 
{k, 1, 4}, {l, 1, 4}], {i, 1, 4}, {j, 1, 4}], Reals]; F2="F"^μσ -> MatrixForm[f]; Grid[{{F1,"  ", F2}}]

"Christoffel Symbole"
christoffel=Simplify[Table[(1/2)Sum[(inversemetrik[[i, s]])
(D[metrik[[s, j]], x[[k]]]+D[metrik[[s, k]], x[[j]]] -D[metrik[[j, k]], x[[s]]]), {s, 1, n}], 
{i, 1, n}, {j, 1, n}, {k, 1, n}]];
christoffelsymbole=Table[If[UnsameQ[christoffel[[i, j, k]], 0], 
{ToString[Γ[i, j, k]], christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1, j}];
TableForm[Partition[DeleteCases[Flatten[christoffelsymbole], Null], 2]]

"Riemann Tensor"
riemann=FullSimplify[Table[
D[christoffel[[i, j, l]], x[[k]]] - D[christoffel[[i, j, k]], x[[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}]];
R=Table[If[UnsameQ[riemann[[i, j, k, l]], 0],
{ToString[R[i, j, k, l]], riemann[[i, j, k, l]]}],
{i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, k - 1}];
TableForm[Partition[DeleteCases[Flatten[R], Null], 2]]

"Ricci Tensor"
ricci=FullSimplify[Table[
Sum[riemann[[i, j, i, l]], {i, 1, n}], {j, 1, n}, {l, 1, n}]];
Ř1=Subscript["Ř", μσ] -> MatrixForm[ricci];
Ricci=FullSimplify[Table[Sum[
inversemetrik[[i, k]] inversemetrik[[j, l]] ricci[[k, l]], {k, 1, 4}, {l, 1, 4}], 
{i, 1, 4}, {j, 1, 4}], Reals];
Ř2="Ř"^μσ -> MatrixForm[Ricci];
Grid[{{Ř1,"  ", Ř2}}]

"Ricci Skalar"
Ř=Simplify[Sum[inversemetrik[[i, j]] ricci[[i, j]], {i, 1, n}, {j, 1, n}]]; "Ř"->Ř

rplc[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]->
t'[τ])/.Derivative[1][r[τ]]->r'[τ])/.Derivative[1][θ[τ]]->θ'[τ])/.Derivative[1][φ[τ]]->φ'[τ]; rple[y_]:=(((((((y/.t->t[τ])/.r->r[τ])/.θ->θ[τ])/.φ->φ[τ])/.Derivative[1][t[τ]]-> 
t'[τ])/. r´->r'[τ])/. θ´->θ'[τ])/. φ´->φ'[τ];list[y_]:=y[[1]]==y[[2]];

"Zeitdilatation"
ṫ=rple[Simplify[Normal[Solve[
-μ==gtt (ť^2)+grr ((r´)^2)+gθθ ((θ´)^2)+gφφ ((φ´)^2) + 2 gtφ ť φ´, ť, Reals]]][[2, 1, 2]]]-щ;
"t'[τ]"==ṫ

"Bewegungsgleichungen"
geodäsie=Simplify[Table[-Sum[
christoffel[[i, j, k]] x[[j]]' x[[k]]'+q f[[i, k]] x[[j]]' metrik[[j, k]], 
{j, 1, n}, {k, 1, n}], {i, 1, n}]];bewegungsgleichung=Table[{x[[i]]''[τ]==Simplify[rplc[geodäsie[[i]]], Reals]}, {i, 1, n}];

geodesic1=bewegungsgleichung[[1]][[1]]
geodesic2=bewegungsgleichung[[2]][[1]]
geodesic3=bewegungsgleichung[[3]][[1]]
geodesic4=bewegungsgleichung[[4]][[1]]