Needs["VariationalMethods`"]

(*Find Lagrangian by computing the metric from the line element*)
Lagrange = Module[{}, z = f[x[t], y[t]]; r = {x[t], y[t], z};
   ds = Sqrt[D[r, t] . D[r, t]];
   ds*ds
   ];
Print["Lagrange=", Lagrange]
(*Mathematica differentiates everything and writes down the Euler \
Lagrange Equation for us*)
{eq1, eq2} = EulerEquations[Lagrange, {x[t], y[t]}, {t}];

(*Solve it,nb:since we only care about the distance we can put any \
time interval*)

S = NDSolve[{eq1, eq2, x[0] == x0, x[1] == x1, y[0] == y0, 
    y[1] == y1}, {x[t], y[t]}, {t, 0, 1}];
xsol[t_] := Evaluate[x[t] /. S[[1, 1]]]
ysol[t_] := Evaluate[y[t] /. S[[1, 2]]]
zsol[t_] := Evaluate[f[xsol[t], ysol[t]]];
solution = 
  ParametricPlot3D[{xsol[t], ysol[t], f[xsol[t], ysol[t]]}, {t, 0, 1},
    PlotRange -> All, PlotStyle -> Red];
Show[surface, solution, ptsPlot]
distance = 
 NIntegrate[Sqrt[xsol'[t]^2 + ysol'[t]^2 + zsol'[t]^2], {t, 0, 1}] // 
  NumberForm[#, 10] &(*Arc Length*)

Second Edit:

I found examples where this notebook fails to give me the shortest distance, for example use:

(* Define Surface*)
f[x_, y_] := 1/(x^2 + y^2 + 0.2)
surface = Plot3D[f[x, y], {x, -1, 1}, {y, -1, 1}];

(*Determine points*)
{x0, y0} = {-1, -1};
{x1, y1} = {1, 1};
ptsPlot = 
  ListPointPlot3D[{{x0, y0, f[x0, y0]}, {x1, y1, f[x1, y1]}}, 
   PlotRange -> All, PlotStyle -> {Black, PointSize[0.05]}];
Show[surface, ptsPlot]

This tells me that the trajectory is by climbing the summit of the mountain, which must be clearly wrong. Is this an error in the NDSolve?