The upshot of this answer is as follows: if a path satisfies the Euler-Lagrange equations for $L^2/2$, then it will satisfy the Euler-Lagrange equations for $L$, but the converse does not hold unless the path has affine parameterization.
Let $L = L(x, \dot x)$ be a lagrangian that is a local function of only position and velocity, then a parameterized path $x(s) = (x^i(s))$ on $M$ is said to satisfy the Euler-Lagrange equations for $L$ provided
\begin{align}
\frac{\partial L}{\partial x^i}(x(s), \dot x(s)) - \frac{d}{ds}\frac{\partial L}{\partial \dot x^i}(x(s), \dot x(s)) = 0
\end{align}
for all $i$ and for all $s$ in the domain of $x$.
Lemma 1. If $x$ satisfies the Euler-Lagrange equations for $L$, then the Beltrami Identity holds for $x$:
$$
\frac{d}{ds}L(x(s), \dot x(s)) = \frac{d}{ds}\left(\frac{\partial L}{\partial \dot x^i}\big(x(s), \dot x(s)\big)\cdot \dot x^i(s)\right)
$$
for all $s$ in the domain of $x$.
Proof. Try it yourself! The proof hinges on the fact that $L$ is a local function of only $x$ and $\dot x$.
Lemma 2. If $L(x,\dot x) = \sqrt{g_{ij}(x)\dot x^i\dot x^j}$, then $L$ satisfies the following identity:
$$
\frac{\partial (L^2/2)}{\partial \dot x^i}(x, \dot x) \dot x^i = L(x,\dot x)^2
$$
Proof. Try this yourself too!
Corollary. If $L(x,\dot x) = \sqrt{g_{ij}(x)\dot x^i\dot x^j}$, and $x$ satisfies the Euler-Lagrange equations for $L^2/2$, then $x$ satisfies the Euler-Lagrange equations for $L$.
Proof. If $x$ satisfies the Euler-Lagrange equations for $L^2$, then Lemma 1 gives the following Beltrami identity (we use notational shorthand here -- all expressions should be evaluated on $x(s)$)
$$
\frac{d(L^2/2)}{ds} = \frac{d}{ds} \frac{\partial (L^2/2)}{\partial \dot x^i}\cdot \dot x^i
$$
On the other hand, evaluating both sides of Lemma 2 on $x(s)$, and taking the derivative of both sides with respect to $s$ gives
$$
\frac{d}{ds} \frac{\partial (L^2/2)}{\partial \dot x^i}\cdot \dot x^i = \frac{d(L^2)}{ds}
$$
Combining these facts shows that $d(L^2)/ds = 0$ which implies that $L^2$ is constant along $x(s)$ and therefore that $L$ is also constant along $x(s)$:
$$
\frac{dL}{ds} = 0.
$$
Now, we separately notice that since $x$ satisfies the Euler-Lagrange equations for $L^2/2$, we have
\begin{align}
0
&= \frac{\partial(L^2/2)}{\partial x^i} - \frac{d}{ds} \frac{\partial (L^2/2)}{\partial \dot x^i} \\
&= L\left(\frac{\partial L}{\partial x^i} - \frac{d}{ds}\frac{\partial L}{\partial \dot x^i}\right) - \frac{dL}{ds}\frac{\partial L}{\partial \dot x^i} \tag{$\star$}\\
&= L\left(\frac{\partial L}{\partial x^i} - \frac{d}{ds}\frac{\partial L}{\partial \dot x^i}\right)
\end{align}
and therefore as long as $L\neq 0$, we see that $x$ satisfies the Euler-Lagrange equations for $L$ as was desired.
The crucial point here is that because of the specific form of $L$, any path satisfying the Euler-Lagrange equation for $L^2/2$ has the nice property that $dL/ds = 0$ along the path. This allows one to kill the term in $(\star)$ which is the term that is the essential difference between the Euler-Lagrange equations for $L^2/2$ and the Euler-Lagrange equations for $L$.
However, if $x$ satisfies the Euler-Lagrange equations for $L$, then it is not necessarily the case that $dL/ds = 0$ along $x$, so in this case, one can't kill that term in $(\star)$, so it need not be a solution to the Euler-Lagrange equation for $L^2/2$.
Nonetheless, if $x$ is affinely parameterized, then it will automatically have the property that $L$ is constant along it, so it will automatically satisfy both Euler-Lagrange equations.
In fact, using parts of the computations above, it is not hard to show that
Proposition. Let $L(x, \dot x) = \sqrt{g_{ij}(x)\dot x^i\dot x^j}$. A path $x$ is an affinely parameterized geodesic if and only if it solves the Euler-Lagrange equations of both $L$ and $L^2/2$.
So the Euler-Lagrange equations of $L^2/2$ yield all affiniely parameterized geodesics, while the Euler-Lagrange equations of $L$ yield all geodesics, regardless of parameterization.
Best Answer
They are all equivalent. The answer to your other question is: the Hamiltonian approach usually works best.
Geodesics can be defined a few ways, since the connection of spacetime is taken to be Levi-Civita.
Let $(\mathcal{M},g)$ denote spacetime $\mathcal{M}$ with metric $g$ and $\gamma:\mathbb{R}\to M,t\mapsto \gamma(t)$ be a curve on $\mathcal{M}$. If $\nabla$ is the Levi-Civita connection on spacetime, then the geodesic equation is $$\nabla_{\dot\gamma}\dot\gamma=0,$$ where $\dot\gamma$ is the tangent vector of $\gamma$. Some manipulations can bring this into the form given in the OP. On the other hand, we can define the length functional $$\ell [\gamma]:=\int_a^b\sqrt{-g(\dot\gamma,\dot\gamma)}\,\mathrm{d}t.$$ (The negative appears because $\gamma$ is usually taken to be timelike for GR purposes.) Then it may be shown that $$\frac{\delta\ell}{\delta\gamma}=0\Longleftrightarrow\nabla_{\dot\gamma}\dot\gamma=0.$$ So the geodesic problem of GR is actually an Euler-Lagrange type problem with Lagrangian $L=\sqrt{-g(\dot\gamma,\dot\gamma)}$.
We can also get a Hamiltonian approach. We first note that the energy functional $$E[\gamma]=\frac{1}{2}\int_a^bg(\dot\gamma,\dot\gamma)\,\mathrm{d}t$$ has identical Euler-Lagrange equations as the length functional. Proof: Let $D$ be the operator $$Df=\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial f}{\partial \dot x}-\frac{\partial f}{\partial x}$$ where $f=f(x,\dot x,t)$, so that $Df=0$ is the Euler-Lagrange equation for $f$. Then a short computation shows that $D\ell=\ell\cdot DE$, so $D\ell=0\implies DE=0$. We note that $H:=\frac{1}{2}g(\dot\gamma,\dot\gamma)$ is reminiscent of the $\frac{1}{2}m\mathbf{v}^2$ term of the classical free particle Hamiltonian. (Hence the name energy functional.)
We can do this more fancily on the cotangent bundle $T^*\mathcal{M}$. Locally trivialize the cotangent bundle in a chart, so we have coordinates $(x^\mu,p_\mu)$. Then put $$H(x,p):=\frac{1}{2}g^{\mu\nu}(x)p_\mu p_\nu.$$ The Hamiltonian equations $$\dot x^\mu=\frac{\partial H}{\partial p_\mu},\quad \dot p^\mu=-\frac{\partial H}{\partial x^\mu}$$ are equivalent to the geodesic equation. The flow obtained from these equations is a Hamiltonian flow.
It turns out that the first Hamilton approach is very useful. It is quite easy to vary the integral $E[\gamma]$. This, combined with the Killing method for obtaining first integrals of the geodesic equation, make for a better strategy than computing the Christoffel symbols and solving the equations outright.
It should be noted that there is another method for solving the geodesic equations: apply the methods of Hamilton-Jacobi theory to the Hamiltonian system described above. This method is used in N. Straumann, General Relativity (2013) to find the geodesics of Kerr spacetime. See also V.I. Arnold, Mathematical Methods of Classial Mechanics (1989) for more on Hamiltonian systems in general.
It should be noted that the Lagrangian and Hamiltonian described above are not actually Legendre transforms of each other.