ClearAll["Global`*"];
Remove["Global`*"];

X[u_, v_]:={u, v, ((u^2)-(v^2))/2}
X[pair_]:=X[pair[[1]], pair[[2]]]
uu[i_]:=If[i==1, u, v];

uut[i_,t_]:= Evaluate[If[i==1, uu1[t], uu2[t]]];

X1[u0_,v0_]:=D[X[u,v], u]/.{u->u0, v->v0}
X2[u0_,v0_]:=D[X[u,v], v]/.{u->u0, v->v0}

dXdui [i_, u0_, v0_]:=If[i==1, X1[u0, v0], If[i==2, X2[u0, v0], 0]];

Nvec[u_, v_]:=Normalize[Cross[X1[u,v], X2[u,v]]];

v1[u_,v_]:=Normalize[X1[u,v]];

v2[u_,v_]:=Cross[Nvec[u,v],v1[u,v]];

DirVec[u_,v_, θ_]:= v1[u,v]Cos[θ]+v2[u,v]Sin[θ];

BinVec[u_,v_,θ_]:=DirVec[u,v, θ + π/2];

NIntPlane[u_, v_, θ_]:=Cross[Normalize[DirVec[u, v, θ]], Nvec[u, v]];

g[i_, j_,u0_, v0_]:= dXdui[i, u0, v0].dXdui[j, u0, v0];
GMat[u0_, v0_]:= {{g[1, 1,u0, v0],g[1,2,u0, v0]},{g[2, 1,u0, v0], g[2, 2,u0, v0]}};
GMatInv[u0_, v0_]:= Inverse[GMat[u0,v0]];

Christoff1st[i_, j_,k_, u0_, v0_]:=(1/2)(D[g[i, k, u0,v0], uu[j]]+D[g[j, k, u0,v0], uu[i]]-D[g[i, j, u0,v0], uu[k]])
Christoff2nd[r_, i_, j_, u0_, v0_]:=Sum[GMatInv[u0,v0][[k, r]] Christoff1st[i, j,k, u0, v0], {k, 1, 2}]

GeodesicEq[r_,t_]:=D[uut[r,t], {t, 2}]+Sum[Sum[Christoff2nd[r, i, j, u, v] D[uut[i,t],t]D[uut[j,t],t], {i, 1, 2}] , {j, 1, 2}]/.{u->uu1[t], v->uu2[t]}//Simplify

(*ics[x0_, y0_, θ_]:={uu1[10]==x0, uu2[10]==y0, uu1'[10]==Cos[θ], uu2'[10]==Sin[θ]};*)
ics[x0_, y0_, θ_]:={uu1[10]==x0, uu2[10]==y0, uu1'[10]==DirVec[x0, y0, θ][[1]], uu2'[10]==DirVec[x0, y0, θ][[2]]};

SystemofEqns[x0_, y0_, θ_]:=Append[Table[GeodesicEq[r, t]==0, {r, 1, 2}],ics[x0, y0, θ]]

w1[x0_, y0_,θ0_]:= Evaluate[DirVec[x0, y0,θ0]];
w2[x0_, y0_]:= Evaluate[Nvec[x0,y0]];

Plane[u_, v_, x0_, y0_, θ_]:=X[x0, y0]+ u w1[x0, y0,θ]+v w2[x0, y0];

bb[t_, x_, y_, θ2_]:=b/.NSolve[X[x0, y0]+t0*w1[x0, y0, θ] + b*w2[x0, y0]==X[x0+(t0*w1[x0, y0, θ] + b*w2[x0, y0])[[1]], y0+(t0*w1[x0, y0, θ] + b*w2[x0, y0])[[2]]],b][[1]]/.{x0->x, y0->y, θ->θ2, t0->t}

α[t_, x0_, y0_, θ_]:=X[x0, y0]+t*w1[x0, y0, θ] + bb[t, x0, y0, θ] *w2[x0, y0]

K:=Abs[((bInterp''[0])/(((1+bInterp'[0]^2))^(3/2)))];

OscCenter := αInterp[0] + If[Re[αInterp''[0].NvecOld] ≥ 0, NvecOld/K, -NvecOld/K];

OscCircle[t_]:=OscCenter+Cos[t]NvecOld/K +Sin[t] DirVecOld/K;

ReinterpolateCurve[x0_, y0_, θ_]:=
{
pts=Table[{t,bb[t, x0, y0, θ]}, {t, 0, 20, .2}];
bInterp=Interpolation[pts];
αInterp[t_] = X[x0, y0]+t*w1[x0, y0, θ] + bInterp[t] *w2[x0, y0];
NvecOld = Nvec[x0, y0];
DirVecOld = DirVec[x0, y0, θ];
θOld=θ;
OscCenterOld = Evaluate[OscCenter[t]];
KOld=Evaluate[K[t]];
}

ComputeGeodesic[x0_, y0_, θ_]:=
{
sol=NDSolve[SystemofEqns[x0, y0, θ], {uu1, uu2}, {t, 0, 20}];
}
ComputeGeodesic[1, 2, 0]

αInterp[t_] = X[x0, y0]+t*w1[x0, y0, θ] + bInterp[t] *w2[x0, y0]/.{x0->1, y0->2,θ->0};

Manipulate[
{Show[{
	ParametricPlot3D[{Plane[u, v, x0, y0, θ]}, {u, -4, 4}, {v,-4, 4}, 
					PlotRange->{{-4, 4}, {-4, 4}, {-4, 4}} ,PlotStyle->Opacity[.1],Mesh->None, BoxRatios->Automatic],
	ParametricPlot3D[{X[u,v]}, {u, -4, 4}, {v,-4, 4},
					PlotRange->{{-4, 4}, {-4, 4}, {-4, 4}} ,PlotStyle->Opacity[.4],                  Mesh->None, BoxRatios->Automatic],
	ParametricPlot3D[X[Flatten[Evaluate[{uu1[t],uu2[t]}/.sol]]], {t, 0, 20},PlotStyle->{Directive[Cyan, Thick]}],
	Graphics3D[{PointSize[.02], Blue, Point[X[x0, y0]]}],
	Graphics3D[{Thick, Green, Line[{X[x0, y0],X[x0, y0]+Nvec[x0,y0]}]}],
         Graphics3D[{Thick, Purple,Line[{ X[x0, y0],X[x0, y0]+DirVec[x0, y0, θ]}]}], 
         Graphics3D[{Thick, Cyan, Line[{X[x0, y0], X[x0, y0]+BinVec[x0, y0,θ]}]}]
}]
Button["Compute Geodesic", ComputeGeodesic[x0, y0, θ]]},
{{x0,1}, -4, 4}, {{y0, 2}, -4, 4}, {θ, 0, 2π}]