This Mathematica code was shown on a slide during the Monday, October 21, 2024 lecture I had left out the "Options" part to squeeze it into 10 lines and to be able to brag about brevity. I feel that brevity is a clarity enhancer. I had programmed this already 30 years ago in 1995 in mathematica 2.1 for this course 1995 but the code has become since much shorter. I used to program by defining functions like r[u_,v_] etc which essentially doubles the length of the code. The DSolve routine still needs functions. The code can easily be adapted to solve the geodesic diff equations in arbitrary dimensions. 
(* Solving the Geodesic diff equations in Mathematica, Oliver Knill, 2024    *)
(* First programmed in Mathematica 2.1 in 1995 at Caltech for                *)
(* https://people.math.harvard.edu/~knill/teaching/math109_1995/geometry.pdf *)
(* This 10 line code for solving the geodesics has been shared with the math *)
(* Math 136 class https://people.math.harvard.edu/~knill/teaching/math136    *)

r={(3+Sin[v]) Cos[u],(3+Sin[v]) Sin[u],Cos[v]+Sin[u+v]/2};
ru=D[r,u];rv=D[r,v];drT={ru,rv}; dr=Transpose[drT]; g=drT.dr; gi=Inverse[g]; 
X={u,v}; Y={ud,vd}; d=2; T=100; 
c[i_,j_,k_]:=(D[g[[j,k]],X[[i]]] +D[g[[k,i]],X[[j]]] -D[g[[i,j]],X[[k]]])/2;
Christoffel[i_,j_,k_]:=Sum[gi[[k,l]]*c[i,j,l],{l,d}];
F=Simplify[Table[-Sum[Christoffel[i,j,k] Y[[i]] Y[[j]], {i,d},{j,d}],{k,d}]];
o=NDSolve[{U''[t] == F[[1]] /. u->U[t] /. v->V[t] /. ud->U'[t] /. vd->V'[t],
           V''[t] == F[[2]] /. u->U[t] /. v->V[t] /. ud->U'[t] /. vd->V'[t],
           U[0]==0.3, V[0]==0.5, U'[0]==0.2, V'[0]==1.6}, {U[t],V[t]},{t,0,T}]; 
SetOptions[ParametricPlot3D,{Mesh->False,Boxed->False,Axes->False}];
Show[{ParametricPlot3D[r,{u,0,2Pi},{v,0,2Pi}],
      ParametricPlot3D[(r /. u -> U[t] /. v -> V[t]) /. o,{t,0,T}]}]