If we have a metric $g_{\mu \nu}$, defined in a Riemannian manifold we can write the equation of the geodesic:
$$\frac{d^2x^\mu}{dt^2}+\Gamma^\mu_{\alpha\beta}\frac{dx^\alpha}{dt}\frac{dx^\beta}{dt}$$ in which $\Gamma^\mu_{\alpha\beta}$ are the Christoffel symbols. The geodesic in that metric can be obtained also using the Euler – Lagrange equations:
$$\frac{d}{dt}\frac{\partial L}{\partial \dot q_k}-\frac{\partial L}{\partial q_k}=0$$
Where the $q_k$ are the generalized coordinates. How can I find the Lagrangian $L$ in order to find the same geodesic obtained using the affine connection symbols?
[Math] Geodesic and Euler – Lagrange equation
calculus-of-variationsdifferential-geometry
Related Solutions
I'm not sure what happens for general curves, but I think I can prove the following:
Let $\gamma:[0,1]\rightarrow M$ be any injective curve segement. Then there is a Riemannian metric for which $\gamma$ is a geodesic. If instead $\gamma$ is a simple closed curve and $\gamma'(0) = \gamma'(1)$, the conclusion still holds.
I'm not sure what happens in the other cases.
Here's the idea of the proof in the (slightly harder) second case:
Pick a background Riemannian metric once and for all. The normal bundle of $\gamma$ embeds into $M$ via the exponential map (for a suitably short time). Call the image of this embedding $W$. Choose an open set $V$ with the property that $V\subseteq \overline{V}\subseteq W$ and let $U = M-\overline{V}$. Notice that $W\cup U = M$, so we can find partition of unity $\{\lambda_U,\lambda_W\}$ subordinate to $\{U,W\}$.
Now, the classification of vector bundles over circles is easy: There are precisely 2 of any rank - the trivial bundle of rank $k$ and the Möbius bundle + trivial bundle of rank $k-1$. The point is that both of these have (flat) metrics where the $0$ section ($\gamma$) is a geodesic.
Since $W$ is diffeomorphic to a vector bundle over the circle, we can assume it has a metric $g_W$ for which $\gamma$ is a geodesic. Now, pick any Riemannian metric $g_U$ on $U$. Finally, define the metric $g_M$ on $M$ by $\lambda_W g_W + \lambda_U g_U$. This is a convex sum of metrics, and hence is a metric. Near $\gamma$, $\lambda_U \equiv 0$ and $\lambda_W\equiv 1$, so the metric near $\gamma$ looks just like $g_W$, so $\gamma$ is a geodesic in $M$.
Implicit coordinates
If everything is given in impicit coordinates, the easiest way to check if a curve is geodesic, is to check if the curvature vector of the curve is perpendicular to the manifold: $$ \pmb k = \ddot{\pmb\gamma}-\dot{\pmb\gamma}\frac{\ddot{\pmb\gamma}\cdot\dot{\pmb\gamma}}{\dot{\pmb\gamma}^2}, \qquad |\pmb k\cdot\pmb n|=|\pmb k||\pmb n|, $$ where $\mathbf n$ is a normal to the manifold.
So for your example: $$ \dot{\pmb\gamma}=(-\sin s, \cos s, 0),\qquad \ddot{\pmb\gamma}=(-\cos s, -\sin s, 0),\\ \pmb k=\ddot{\pmb\gamma}=(-\cos s, -\sin s, 0) $$
So you can see that $\pmb k$ is perpendicular to the sphere.
Local coordinates
However, if you want to use the equations you presented, you need to keep in mind that both equations assume local coordinates $x^\mu$. That is if the manifold is embedded is $\mathbb R^3$ with coordinates $\xi^i = (\xi, \eta, \zeta)$, then local coordinates are some functions $x\equiv x^1=x^1(\xi,\eta,\zeta)$, $y\equiv x^2=x^2(\xi,\eta,\zeta)$. And the curve is assumed to be expressed in this local coordinates.
So as an example, let's consider unit sphere. We can select spherical coordinates as local coordinates: $$ \xi = \cos y\cos x,\qquad \eta=\cos y\sin x, \qquad\zeta = \sin y,\\ x = \begin{cases}\arccos\frac{\xi}{\sqrt{\xi^2+\eta^2}},& \eta\ge 0,\\ -\arccos\frac{\xi}{\sqrt{\xi^2+\eta^2}},& \eta< 0,\end{cases}, \qquad y = \arctan\frac{\zeta}{\sqrt{\xi^2+\eta^2}}, $$
Then we can calulate metric tensor of our manifold in local coordinates: $$ g_{\alpha\beta}=\mathbf e_\alpha\cdot\mathbf e_\beta=\frac{\partial\xi^i}{\partial x^\alpha}\frac{\partial\xi^i}{\partial x^\beta}=\begin{pmatrix}\cos^2y& 0\\0 & 1\end{pmatrix}. $$
By inverting the matrix we obtain the inverse metric tensor: $$ g^{\alpha\beta} = \begin{pmatrix}1/\cos^2y& 0\\0 & 1\end{pmatrix}. $$
Christoffel symbols by defintion: $$ \Gamma^\mu_{\alpha\beta} = \frac{\partial \mathbf e_\alpha}{\partial x^\beta}\cdot \mathbf e^\mu= \frac{\partial^2 \xi^i}{\partial x^\alpha\partial x^\beta}g^{\mu\nu}\frac{\partial \xi^i}{\partial x^\nu}. $$
Computation gives us: $$ \Gamma^1_{\alpha\beta} = \begin{pmatrix}0 & -\tan y\\-\tan y& 0\end{pmatrix},\qquad \Gamma^2_{\alpha\beta} = \begin{pmatrix}\sin y\cos y & 0\\0 & 0\end{pmatrix}. $$
Now let's consider a curve on a sphere with equation $x=s$, $y=y_0$ (some parallel with latitude $y_0$). The geodesic equation will give us: $$ \mu=1:\qquad 0 + \begin{pmatrix}1\\ 0\end{pmatrix}^T\begin{pmatrix}0 & -\tan y\\-\tan y& 0\end{pmatrix}\begin{pmatrix}1\\ 0\end{pmatrix} =0\\ \mu=2: \qquad 0 + \begin{pmatrix}1\\ 0\end{pmatrix}^T\begin{pmatrix}\sin y\cos y & 0\\0& 0\end{pmatrix}\begin{pmatrix}1\\ 0\end{pmatrix} = \sin y\cos y = 0 $$
So first equation is always true, while second is true only if $y_0=0$ which is the equator (we don't consider degraded curves $y_0=\pm\pi/2$).
This the most straightforward way to use geodesics. However, it implies a lot of calculations.
Euler-Lagrange equation
You have made a mistake in Euler-Lagrange equation. The Lagrangian $L=\frac12\dot{\pmb\gamma}^2$, and not just $\dot\gamma$. So the equation should read: $$ \frac12\frac{\partial}{\partial x^\mu}\dot{\pmb\gamma}^2 - \frac12\frac{d}{ds}\frac{\partial\dot{\pmb\gamma}^2}{\partial \dot x^\mu}=0 $$
Now how to use it. First, let's notice that we can rewrite the Lagrangian in local coordinates: $L=\dot{\pmb\gamma}^2=g_{\alpha\beta}\dot x^\alpha\dot x^\beta$. It's just the definition of scalar product in curvilinear coordinate system. Then we need to keep in mind that $x^\mu$ and $\dot x^\mu$ are assumed to be independent. So in first term: $$ \frac12\frac{\partial}{\partial x^\mu}\left(g_{\alpha\beta}\dot x^\alpha\dot x^\beta\right) = \frac12\frac{\partial g_{\alpha\beta}}{\partial x^\mu}\dot x^\alpha\dot x^\beta. $$
The inside part of second term: $$ \frac{\partial g_{\alpha\beta}\dot x^\alpha\dot x^\beta}{\partial \dot x^\mu} = g_{\alpha\beta} \frac{\partial\dot x^\alpha}{\partial \dot x^\mu}\dot x^\beta + g_{\alpha\beta} \frac{\partial\dot x^\beta}{\partial\dot x^\mu}\dot x^\alpha = g_{\alpha\beta}\delta_{\alpha\mu}\dot x^\beta+g_{\alpha\beta}\delta_{\beta\mu}\dot x^\alpha = 2g_{\mu\alpha}\dot x^\alpha. $$ In last part we used the symmetry of metric tensor: $g_{\alpha\beta}=g_{\beta\alpha}$.
So now we can rewrite our equation as: $$ \frac12\frac{\partial g_{\alpha\beta}}{\partial x^\mu}\dot\gamma^\alpha\dot\gamma^\beta -\frac{d}{ds}\left(g_{\mu\alpha}\dot\gamma^\alpha\right). $$
Let's consider the same metric and curve as in previous section. $$ \frac{\partial g_{\alpha\beta}}{\partial x^1} = 0,\qquad \frac{\partial g_{\alpha\beta}}{\partial x^2} = \begin{pmatrix} -2\cos x\sin x&0\\0&0 \end{pmatrix} $$
Covector $g_{\mu\alpha}\dot\gamma^\alpha=\begin{pmatrix}\cos^2 y\\0\end{pmatrix}$ doesn't depend on $s$, so its full derivative $d/ds$ is zero. Thus we end with the same equation: $$\frac12\begin{pmatrix}1\\ 0\end{pmatrix}^T\begin{pmatrix}-2\sin y\cos y & 0\\0& 0\end{pmatrix}\begin{pmatrix}1\\ 0\end{pmatrix}=0 $$
Finally, it's possible to show that the equation is the same: $$ \frac{d}{ds}(g_{\mu\alpha}\dot\gamma^\alpha) = \frac{\partial g_{\mu\alpha}}{\partial x^\beta}\dot\gamma^\beta\dot\gamma^\alpha + g_{\mu\alpha}\ddot\gamma^\alpha = \frac12\left(\frac{\partial g_{\mu\alpha}}{\partial x^\beta}+\frac{\partial g_{\mu\beta}}{\partial x^\alpha}\right)\dot\gamma^\beta\dot\gamma^\alpha + g_{\mu\alpha}\ddot\gamma^\alpha. $$ In last part we used the symmetry of $\dot\gamma^\beta\dot\gamma^\alpha$ and split the term in two and renamed indices in one of the new terms. Thus, we can write our equation as: $$ g_{\mu\alpha}\ddot\gamma^\alpha +\frac12\left(\frac{\partial g_{\mu\alpha}}{\partial x^\beta}+\frac{\partial g_{\mu\beta}}{\partial x^\alpha} - \frac{\partial g_{\alpha\beta}}{\partial x^\mu}\right)\dot\gamma^\beta\dot\gamma^\alpha = g_{\mu\alpha}\ddot\gamma^\alpha + \Gamma_{\mu\,\alpha\beta}\dot\gamma^\beta\dot\gamma^\alpha = 0 $$ Finally, we multiply the equation by $g^{\mu\nu}$ and get: $$ \ddot\gamma^\nu + \Gamma^\nu_{\alpha\beta}\dot\gamma^\beta\dot\gamma^\alpha = 0 $$
Best Answer
The Lagrangian is the standard measurement of a distance. That is $$ L=-\kappa\int ds $$ being $\kappa>0$ a constant, that in your case is $$ S=-\kappa\int\sqrt{g_{\mu\nu}dx^\mu dx^\nu}. $$ Then you are supposed to parametrize the coordinates as $x_\mu=x_\mu(t)$ and the computation becomes a standard one. Notice that, for a Minkowski metric, this reduces to the action of a relativistic free particle.