Curvature and the geodesic equations

Copy/paste the following cell in a notebook and edit the data between the red lines. Then evaluate it to see the results.
Details on the code and how it works can be found at the end of the notebook.

Clear[coord, metric,inversemetric, affine, riemann, lriemann, uriemann, ricci, scalar, einstein,weyl, geodesic,R,G,τ, i,j,k,l,s];
Clear[r,θ,ϕ,t,χ,a,m,M]; 

(*----------------------------------------------------------------------------*)
(* This is what you need to set: *)
coord = {v,r,θ,ϕ};
n          =Length[coord];
metric={
{-(1-2(M/r)),1                ,0                         ,                                        0},
{1                  ,0                ,0                         ,                                        0},
{0                 ,0                 ,(r^2)                        ,                                         0},
{0                 ,0                 ,0                          ,(r^2) ((Sin[θ])^2)                      }};
(*----------------------------------------------------------------------------*)

inversemetric=FullSimplify[Inverse[metric]];
Print["-----------------------------------------------------------------------"];
Print["The Manifold has dimension   n= ",n,"\nCoordinate system:          ",coord];
Print["-----------------------------------------------------------------------"];
Print["(Subscript[g, μν])=",metric//MatrixForm];
Print["(g^(μν))=",inversemetric//MatrixForm];
Print["g  =",Det[metric]//FullSimplify];
affine:=affine=FullSimplify[Table[
(1/2)*Sum[
(*              (g^is) ((∂_k)(Subscript[g, sj])+(∂_j)(Subscript[g, sk])-(∂_s)(Subscript[g, jk]))           *)
(inversemetric[[i,s]])*
(D[metric[[s,j]],coord[[k]] ]+D[metric[[s,k]],coord[[j]] ]-D[metric[[j,k]],coord[[s]] ]),
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n}] ];
(*The non zero Christoffel symbols are computed and selected below: *)
listaffine:=Table[
If[
UnsameQ[affine[[i,j,k]],0],
{Subscript[Superscript[Γ,i-1],j-1,k-1],affine[[i,j,k]]}
  ] ,
{i,1,n},{j,1,n},{k,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Christoffel Symbols:"];
Print[TableForm[Partition[DeleteCases[Flatten[listaffine],Null],2],TableSpacing->{2,2}]];
riemann:=riemann=FullSimplify[Table[
(* (Subscript[(R^i), jkl])=    (∂_k)(Subscript[(Γ^i), lj])                   -         (∂_l)(Subscript[(Γ^i), kj])*)
D[       affine[[i,l,j]],  coord[[k]]           ]-D[affine[[i,k,j]], coord[[l]]          ]+
(*         (Subscript[(Γ^i), ks])             (Subscript[(Γ^s), lj])     -         (Subscript[(Γ^i), ls])            (Subscript[(Γ^s), kj])      *)
Sum[affine[[i,k,s]]affine[[s,l,j]] -    affine[[i,l,s]] affine[[s,k,j]] ,
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}] ];
listriemann:=Table[
If[
UnsameQ[riemann[[i,j,k,l]],0],
{Subscript[Superscript[R,i-1],j-1,k-1,l-1],riemann[[i,j,k,l]]}
] ,
{i,1,n},{j,1,n},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listriemann],Null],2],TableSpacing->{2,2}]];
lriemann:=lriemann=FullSimplify[Table[
Sum[metric[[i,ii]] riemann[[ii,j,k,l]],{ii,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listlriemann:=Table[
If[
UnsameQ[lriemann[[i,j,k,l]],0],
{Subscript[R,i-1,j-1,k-1,l-1],lriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Contravariant Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listlriemann],Null],2],TableSpacing->{2,2}]];
uriemann:=uriemann=FullSimplify[Table[
Sum[
inversemetric[[j,jj]] inversemetric[[k,kk]] inversemetric[[l,ll]] riemann[[i,jj,kk,ll]],{jj,1,n},{kk,1,n},{ll,1,n}
],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listuriemann:=Table[
If[
UnsameQ[uriemann[[i,j,k,l]],0],{Superscript[Superscript[Superscript[Superscript[R,i-1],j-1],k-1],l-1],uriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Covariant Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listuriemann],Null],2],TableSpacing->{2,2}]];
r2=FullSimplify[Sum[lriemann[[i,j,k,l]]uriemann[[i,j,k,l]],{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
Print["-----------------------------------------------------------------------"];
Print["(R^2)= ",r2];
ricci:=ricci=FullSimplify[Table[
Sum[
riemann[[i,j,i,l]],
{i,1,n}
],{j,1,n},{l,1,n}] ];
listricci:=Table[
If[
UnsameQ[ricci[[j,l]],0],
{Subscript[R,j-1,l-1],ricci[[j,l]]}
] ,{j,1,n},{l,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Ricci Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listricci],Null],2],TableSpacing->{2,2}]];
scalar=FullSimplify[Sum[inversemetric[[i,j]]ricci[[i,j]],{i,1,n},{j,1,n}] ];
Print["-----------------------------------------------------------------------"];
Print["Curvature Scalar:"];
Print["R= ",scalar];
einstein:=einstein=FullSimplify[ricci-(1/2)scalar*metric];
listeinstein:=Table[
If[
UnsameQ[einstein[[j,l]],0],
{Subscript[G,j,l],einstein[[j,l]]}
] ,{j,1,n},{l,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Einstein Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listeinstein],Null],2],TableSpacing->{2,2}]];
weyl:=weyl=FullSimplify[Table[
If[n>3,
lriemann[[i,j,k,l]]
-         (1/(n-2))           (metric[[i,k]] ricci[[l,j]] -metric[[i,l]] ricci[[k,j]]-metric[[j,k]] ricci[[l,i]]+metric[[j,l]] ricci[[k,i]])
+(1/((n-1)(n-2)))(metric[[i,k]] metric[[l,j]]-metric[[i,l]] metric[[k,j]]) scalar
(*else, if n<= 3 return 0:*),0],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listweyl:=Table[
If[
UnsameQ[weyl[[i,j,k,l]],0],
{Subscript[ℭ,i-1,j-1,k-1,l-1],weyl[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Weyl Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listweyl],Null],2],TableSpacing->{2,2}]];
geodesic:=geodesic=Simplify[Table[-Sum[affine[[i,j,k]]u[j]u[k],{j,1,n},{k,1,n}],{i,1,n}]];
subst = Table[u[i]-> Subscript[coord[[i]],τ],{i,1,n}];
nlistgeodesic:=Table[{Subscript[coord[[i]],ττ],"+",-geodesic[[i]]/.subst,"= 0"},{i,1,n}];
Print["-----------------------------------------------------------------------"];
Print["Geodesic Equations:"];
Print[TableForm[nlistgeodesic,TableSpacing->{2}]];

Postscript

The code in this notebook is based on James Hartle’s Mathematica programs written as companions to his book “Gravity: An Introduction to Einstein’s General Relativity”, available from the Book’s site:
http://web.physics.ucsb.edu/~gravitybook/mathematica.html
The original program was written by Leonard Parker, University of Wisconsin, Milwaukee (see the Acknowledgement section in the end of this notebook). The code in this notebook is a modification of the original program, adapted for the needs of the course “General Relativity and Cosmology” offered at the National Technical University of Athens, by Konstantinos Anagnostopoulos ( http://physics.ntua.gr/konstant ).The site of the course for the Spring 2023 can be found at  http://physics.ntua.gr/konstant/GR .

This is the Mathematica notebook Curvature and the Einstein Equation available from the book website.  From a given metric (Subscript[g, αβ]) , it computes the components of the following: the inverse metric, (g^(λσ)), the Christoffel symbols or affine connection,

(Subscript[(Γ^λ), μν])=(1/2)(g^(λσ))((∂_μ) (Subscript[g, σν])+(∂_ν) (Subscript[g, σμ])-(∂_σ) (Subscript[g, μν])),

( (∂_(α  ))stands for the partial derivative ∂/∂(x^α)), the Riemann tensor,

(Subscript[(R^λ), μνσ])=(∂_ν) (Subscript[(Γ^λ), μσ])-(∂_σ) (Subscript[(Γ^λ), μν])+(Subscript[(Γ^λ), ην])(Subscript[(Γ^η), μσ])-(Subscript[(Γ^λ), ησ])(Subscript[(Γ^η), μν]),

the Ricci tensor

(Subscript[R, μν])=(Subscript[(R^λ), μλν]),

the scalar curvature,

R=(g^(μν))(Subscript[R, μν]),

and the Einstein tensor,

(Subscript[G, μν])=(Subscript[R, μν])-(1/2)(Subscript[g, μν])R.

You must input the covariant components of the metric tensor (Subscript[g, μν]) by editing the relevant input line in this Mathematica notebook. You may also wish to change the names of the coordinates. Only the nonzero components of the above quantities are displayed as the output. All the components computed are in the coordinate basis in which the metric was specified.