Deriving the Schwarzschild solution
Last edited by Jens Boos, Feb 14, 2024.
Based on this Mathematica script idea: https://mathematica.stackexchange.com/questions/8895/how-to-calculate-scalar-curvature-ricci-tensor-and-christoffel-symbols-in-mathe

(* Schwarzschild metric *)
Clear[g,xx,f,ds2]
xx={t,r,θ,ϕ};
ds2=-f[r] dt^2+dr^2/f[r]+r^2(dθ^2+Sin[θ]^2dϕ^2);
gtt=(1/2)D[ds2,dt,dt]//Simplify;
grr=(1/2)D[ds2,dr,dr]//Simplify;
gθθ=(1/2)D[ds2,dθ,dθ]//Simplify;
gϕϕ=(1/2)D[ds2,dϕ,dϕ]//Simplify;
g={{gtt,0,0,0},{0,grr,0,0},{0,0,gθθ,0},{0,0,0,gϕϕ}};
gInv=Inverse[g];

(* all kinds of geometrical quantities *)
InverseMetric[g_]:=Simplify[Inverse[g]];
ChristoffelSymbol[g_,xx_]:=Block[{n,ig,res},n=4;ig=InverseMetric[g];
res=Table[(1/2)*Sum[ig[[i,s]]*(-D[g[[j,k]],xx[[s]]]+D[g[[j,s]],xx[[k]]]+D[g[[s,k]],xx[[j]]]),{s,1,n}],{i,1,n},{j,1,n},{k,1,n}];
Simplify[res]]
RiemannTensor[g_,xx_]:=Block[{n,Chr,res},n=4;Chr=ChristoffelSymbol[g,xx];
res=Table[D[Chr[[i,k,m]],xx[[l]]]-D[Chr[[i,k,l]],xx[[m]]]+Sum[Chr[[i,s,l]]*Chr[[s,k,m]],{s,1,n}]-Sum[Chr[[i,s,m]]*Chr[[s,k,l]],{s,1,n}],{i,1,n},{k,1,n},{l,1,n},{m,1,n}];
Simplify[res]]
RicciTensor[g_,xx_]:=Block[{Rie,res,n},n=4;Rie=RiemannTensor[g,xx];
res=Table[Sum[Rie[[s,i,s,j]],{s,1,n}],{i,1,n},{j,1,n}];
Simplify[res]]
RicciScalar[g_,xx_]:=Block[{Ricc,ig,res,n},n=4;Ricc=RicciTensor[g,xx];ig=InverseMetric[g];
res=Sum[ig[[s,i]] Ricc[[s,i]],{s,1,n},{i,1,n}];
Simplify[res]]
EinsteinTensor[g_,xx_]:=Block[{RicT,RicS,res,n},n=4;RicT=RicciTensor[g,xx];RicS=RicciScalar[g,xx];
res=Table[RicT[[i,j]]-1/2RicS g[[i,j]],{i,1,n},{j,1,n}];
Simplify[res]
]
WeylTensor[g_,xx_]:=Block[{Rie,Ric,R,ig,res,n},n=4;Rie=RiemannTensor[g,xx];Ric=RicciTensor[g,xx];
R=RicciScalar[g,xx];
res=Table[Sum[g[[i,a]]Rie[[a,j,k,l]],{a,1,n}]-1/2g[[i,k]]Ric[[j,l]]+1/2g[[i,l]]Ric[[j,k]]-1/2g[[j,l]]Ric[[i,k]]+1/2g[[j,k]]Ric[[i,l]]+1/6 R (g[[i,k]]g[[j,l]]-g[[i,l]]g[[j,k]]),{i,1,n},{j,1,n},{k,1,n},{l,1,n}];
Simplify[res]]
RiemannSquared[g_,xx_]:=Block[{Riem,ig,n,res},n=4;Riem=RiemannTensor[g,xx];ig=InverseMetric[g];
res=Sum[g[[i,p]]ig[[j,q]]ig[[k,u]]ig[[l,v]]Riem[[i,j,k,l]]Riem[[p,q,u,v]],{i,1,n},{j,1,n},{k,1,n},{l,1,n},{p,1,n},{q,1,n},{u,1,n},{v,1,n}];
Simplify[res]]
WeylSquared[g_,xx_]:=Block[{Weyl,ig,n,res},n=4;Weyl=WeylTensor[g,xx];ig=InverseMetric[g];
res=Sum[ig[[i,p]]ig[[j,q]]ig[[k,u]]ig[[l,v]]Weyl[[i,j,k,l]]Weyl[[p,q,u,v]],{i,1,n},{j,1,n},{k,1,n},{l,1,n},{p,1,n},{q,1,n},{u,1,n},{v,1,n}];
Simplify[res]]
RicciSquared[g_,xx_]:=Block[{Ric,ig,n,res},n=4;Ric=RicciTensor[g,xx];ig=InverseMetric[g];
res=Sum[ig[[i,k]]ig[[j,l]]Ric[[i,j]]Ric[[k,l]],{i,1,n},{j,1,n},{k,1,n},{l,1,n}];
Simplify[res]]
EinsteinEQ[g_,xx_]:=Block[{RicT,RicS,res,n},n=4;RicT=RicciTensor[g,xx];RicS=RicciScalar[g,xx];
res=Table[RicT[[i,j]]-1/2RicS g[[i,j]]+Λ g[[i,j]],{i,1,n},{j,1,n}];
Simplify[res]
]

(* let's solve the vacuum field equations *)
Gμν=EinsteinTensor[g,xx]

sol=DSolve[Gμν[[3]][[3]]==0,f[r],r]

(* verify *)
f[r_]:=1-2GM/r;
Gμν

(* Kretschmann scalar *)
WeylSquared[g,xx]