(* Note that this is the Schwarzschild metric and not the Weak field approximation *)

g=(({{-(1-2 ((G M)/r)), 0, 0, 0}, {0, (((1-2 ((G M)/r)))^(-1)), 0, 0}, {0, 0, (r^2), 0}, {0, 0, 0, (r^2)((Sin[θ])^2)}}));
ig=Inverse[g];

v={t,r,θ,ϕ};

Christoffel[o1_,u1_,u2_]:=(1/2)(∑_(i=1)^(Length[v]))ig[[o1]][[i]]((∂_(v[[u1]]))g[[i]][[u2]]+(∂_(v[[u2]]))g[[i]][[u1]]-(∂_(v[[i]]))g[[u1]][[u2]])//Simplify

Do[
Do[
Do[
Print[{m,a,b}," ", Christoffel[m,a,b]]
,
{b,1,a}
]
,
{a,1,Length[v]}]
,{m,1,Length[v]}]

vtau=(#[τ]&/@v)
replaceVbyVtau=Rule[#[[1]],#[[2]]]&/@Transpose[{v,vtau}]

Do[
Print[(∂_({τ,2}))vtau[[m]]+(∑_(i=1)^(Length[v]))(∑_(j=1)^(Length[v]))Christoffel[m,i,j]((∂_τ)vtau[[i]])((∂_τ)vtau[[j]])],
{m,1,Length[v]}]

(* Make Ansatz for circular orbit *)

ansatz={τ,R,π/2,vel τ};
replaceByAnsatz=Rule[#[[1]][τ],#[[2]]]&/@Transpose[{v,ansatz}]

Do[
Print[((∂_({τ,2}))ansatz[[m]]+(∑_(i=1)^(Length[v]))(∑_(j=1)^(Length[v]))(Christoffel[m,i,j]/.replaceVbyVtau)((∂_τ)ansatz[[i]])((∂_τ)ansatz[[j]]))/.replaceByAnsatz],
{m,1,Length[v]}]

vel/.Solve[-((G M (2 G M-R))/(R^3))+(2 G M-R) (vel^2)==0,vel]//Simplify//PowerExpand

(* This is the angular velocity of a particle on a circular orbit with radius R around a mass M *)

(* Make Ansatz for general geodesic in the plane θ=π/2 *)

ansatz={t[τ],r[τ],π/2,ϕ[τ]};
replaceByAnsatz=Rule[#[[1]][τ],#[[2]]]&/@Transpose[{v,ansatz}]

Do[
Print[(((∂_({τ,2}))ansatz[[m]]+(∑_(i=1)^(Length[v]))(∑_(j=1)^(Length[v]))(Christoffel[m,i,j]/.replaceVbyVtau)((∂_τ)ansatz[[i]])((∂_τ)ansatz[[j]]))/.replaceByAnsatz)],
{m,1,Length[v]}]