# How to eliminate the error in the code of calculating geodesic curve

I find the code to calculate the geodesic on a general surface from here.

But there is a mistake in calculating the geodesic between point {0, 2 ,f[0,2]} and point {0, -2 ,f[0,-2]} on surface $$f(x,y)=x^2+y^2$$:

f[{u_, v_}] := {u, v, u^2 + v^2}
Needs["VariationalMethods"]
eq = EulerEquations[Sqrt[Total[D[f[{u, v[u]}], u]^2]], v[u], u];
geodesic[{{u1_, v1_}, {u2_, v2_}}] :=
Module[{start, g, sol},
If[u2 < u1, Return[geodesic[{{u2, v2}, {u1, v1}}]]];
sol = ParametricNDSolve[Flatten[{eq, v[0] == v1, v'[0] == a}],
v, {u, 0, u2 - u1}, {a}];
start = a /. FindRoot[Evaluate[(v[a][u2 - u1] - v2 /. sol)], {a, 0}];
g = v[start] /. sol;
Function[t, {u1 + t*(u2 - u1), g[t*(u2 - u1)]}]]
pts = {{0, 2}, {0, -2}}
Show[ParametricPlot3D[
f[{u, v}], {u, -π, π}, {v, -π, π},
PlotStyle -> White, ImageSize -> 500],
ParametricPlot3D[Evaluate[f[geodesic[pts][t]]], {t, 0, 1},
PlotStyle -> Red], BoxRatios -> {1, 1, 1}]


But pts = {{0, 2}, {2, 2}} will not report errors. What can I do to eliminate these errors and get the right results quickly.

• The problem is that you parameterize the geodesics as graph over the u-variable in the u-v-plane. But the projection of the geodesic from {0, 2, f[0,-2]} to {0, 2, f[0,-2]}` onto the u-v-plane obviously cannot be a graph over the u-variable. – Henrik Schumacher Feb 10 at 8:00