steps for convenience of the reader. Setting the variation of the length functional to zero we get $$L\frac{d}{dr}\left(\frac{1}{L}2g_{ij}\dot x^i \right)+g_{ij,k} \dot x^i \dot x^j=0$$ with $L\equiv \sqrt{g_{ij}(x(r))\dot x^i\dot x^j}$ and where we use Einstein summation. Parameterizing by arclength, i.e. setting $L\equiv 1$ (or $r=s$) we get $$2g_{ij}\ddot x^i +2g_{ij,l}\dot x^l\dot x^i +g_{ij,k} \dot x^i \dot x^j=0$$ which, upon multiplcation by $g^{mj}$ (entries of inverse of $(g_{mj})$) and rearranging becomes the geodesic equation ,. In order to show on the other hand that the squared length functional satisfies we may consider the length functional $$l(r,x,s,y)=\int_r^s L\left( x(u),\dot x(u)\right) du$$ and invoke Hamilton-Jacobi calculus. This is done by introducing the variables $p_i=L_{\dot x^i}$, and the associated Hamiltonian defined by $$H(x,p)=\dot x^ip_i-L(x,\dot x).$$ (here and henceforth we use Einstein summation if convenient). Then we may write $$x(t)\equiv x(t;r,x,s,y) \mbox{ and } p(t)\equiv p(t;r,x,s,y),$$ where $x(r;r,x,s,y)=x$ and $x(s;r,x,s,y)=y$. and compute $$\label{1ls}
l_s=-H(x(s),p(s)).$$ Then we may connect $p$ to $l_{y^k}$ by computing $$\label{ly}
\begin{array}{ll}
l_{y^k}=\int_r^s \left( \frac{\partial \dot x^i}{\partial y^k}p_i+\dot x^i\frac{\partial p_i}{\partial y^k}
-H_{x^i}\frac{\partial x^i}{\partial y^k}-H_{p_i}\frac{\partial p_i}{\partial y^k}\right) dt\\
\\
\int_r^s \stackrel{{\bf \cdot}}{\left( \frac{\partial x_i}{\partial y^k}p_i\right)}dt=\frac{\partial x_i}{\partial y^k}p_i \Big|^s_r=p_k(s;r,x,s,y).
\end{array}$$ by invoking the canonical system of equations. This leads to $$\frac{\partial l}{\partial s}+\sum_{ij}g^{ij}\frac{\partial l}{\partial y_i}\frac{\partial 