Is there a good way to compute Christoffel Symbols

I'd like to expand on Hans Lundmark's answer because the question keeps recurring.

Let $q^1,\ldots,q^n$ denote the generalized coordinates. Introduce additional velocity variables $\dot{q}^1,\ldots,\dot{q}^n$. If the dot accent is already reserved for other things in your theory, use another accent to avoid confusion.

First, treat all $q^i$ and $\dot{q}^i$ as pairwise independent variables and define, using Einstein summation convention, $$L(q^1,\ldots,q^n;\dot{q}^1,\ldots,\dot{q}^n) = \frac{1}{2} g_{ij}(q^1,\ldots,q^n)\,\dot{q}^i \dot{q}^j\tag{1}$$ To avoid confusion later on, I have explicitly indicated the formal dependencies of the expressions for $g_{ij}$ and $L$. Note that $L$ can immediately be written down when given the first fundamental form.

Now consider a twice differentiable curve, which makes the $q^i$ functions of some independent new parameter $\tau$, and set $\dot{q}^i = \frac{\mathrm{d}q^i}{\mathrm{d}\tau}$.

Proposition. On the curve parameterized with $\tau$ we have $$g^{kh}\left(\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{q}^h}\right) - \frac{\partial L}{\partial q^h}\right) = \ddot{q}^k + \Gamma^k_{\ ij} \dot{q}^i \dot{q}^j \tag{2}$$ You might recognize the right-hand side as the expression constrained by geodesics, and you might recognize the expression wrapped around $L$ in the left-hand side as the form of Euler-Lagrange differential equations, which correspond to variational problems with the Lagrangian $L$. Therefore $(2)$ hints at a variational foundation for the geodesics equation. Such interpretations make $(2)$ easier to memorize or cross-link with other knowledge, but all we actually need is the identity $(2)$. The following observations should provide enough details to prove $(2)$.

We can rewrite $(2)$ to $$\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{q}^h}\right) - \frac{\partial L}{\partial q^h} = g_{hk} \ddot{q}^k + \Gamma_{hij} \dot{q}^i \dot{q}^j \tag{3}$$ The idea now is, for every $h\in\{1,\ldots,n\}$, to take the left-hand side of $(3)$, plug in expressions for the metric in $L$, and rewrite the thing so that it matches the format of the right-hand side, where dotted variables occur only in the places shown. Then you can read off the Christoffel symbols of the first kind, $\Gamma_{h**}$, from the coefficients of the velocity products.

To obtain $(2)$ and $\Gamma^k_{\ **}$ is then just a matter of multiplication with the inverse metric coefficients matrix $((g^{kh}))$, or equivalently, taking the right-hand side expressions of $(3)$ obtained for $h\in\{1,\ldots,n\}$ and then, for each $k\in\{1,\ldots,n\}$, finding a linear combination whose only second derivative with respect to $\tau$ is $\ddot{q}^k$, with coefficient $1$. This is the form given by $(2)$, so the coefficients of the velocity products are then the Christoffel symbols of the second kind, $\Gamma^k_{\ **}$.

Remember, while doing the partial derivatives of $L$, treat the $q^i$ and the $\dot{q}^i$ as independent formal variables. You will do that with more concrete symbol meanings and metric expressions, but in this moderately abstract setting, you can already refine $$\begin{align} \frac{\partial L}{\partial\dot{q}^h} &= g_{hj}\dot{q}^j\tag{4} \\\frac{\partial L}{\partial q^h} &= \frac{1}{2}\frac{\partial g_{ij}}{\partial q^h}\dot{q}^i \dot{q}^j\tag{5} \end{align}$$ However, when doing the $\frac{\mathrm{d}}{\mathrm{d}\tau}$ outside of $L$, stick to the curve and apply the chain rule accordingly: $$\frac{\mathrm{d}g_{hj}}{\mathrm{d}\tau} = \frac{\partial g_{hj}}{\partial q^i}\,\dot{q}^i\tag{6}$$ Now $(4)$, $(5)$, $(6)$ and the Levi-Civita formula $$\Gamma_{hij} = \frac{1}{2}\left( \frac{\partial g_{hj}}{\partial q^i} + \frac{\partial g_{ih}}{\partial q^j} - \frac{\partial g_{ij}}{\partial q^h} \right)$$ can be used to prove $(3)$ and thereby $(2)$. But I will focus on how to apply that proposition.

Example: Spherical coordinates with radius $r$, longitude $\phi$, latitude $\theta$, with $\theta=\frac{\pi}{2}$ at the equator. At index positions, I will write coordinate names instead of digits. The first fundamental form is $$\mathrm{d}s^2 = \mathrm{d}r^2 + (r^2\sin^2\theta)\,\mathrm{d}\phi^2 + r^2\,\mathrm{d}\theta^2$$ Accordingly, the Lagrangian $L$ is $$L = \frac{1}{2}\left(\dot{r}^2 + (r^2\sin^2\theta)\,\dot{\phi}^2 + r^2\,\dot{\theta}^2\right)$$ We now treat $r,\phi,\theta,\dot{r},\dot{\phi},\dot{\theta}$ as independent variables and get $$\begin{align} \frac{\partial L}{\partial\dot{r}} &= \dot{r} &\frac{\partial L}{\partial r} &= (r\sin^2\theta)\,\dot{\phi}^2 + r\,\dot{\theta}^2 \\\frac{\partial L}{\partial\dot{\phi}} &= (r^2\sin^2\theta)\,\dot{\phi} &\frac{\partial L}{\partial\phi} &= 0 \\\frac{\partial L}{\partial\dot{\theta}} &= r^2\,\dot{\theta} &\frac{\partial L}{\partial\theta} &= (r^2\sin\theta\cos\theta)\,\dot{\phi}^2 \end{align}$$ Now we give up the independence, consider some curve paramterized by $\tau$ and obtain $$\begin{align} \frac{\mathrm{d}}{\mathrm{d}\tau} \frac{\partial L}{\partial\dot{r}} &= \ddot{r} \\\frac{\mathrm{d}}{\mathrm{d}\tau} \frac{\partial L}{\partial\dot{\phi}} &= (r^2\sin^2\theta)\,\ddot{\phi} + 2\,(r\sin^2\theta)\,\dot{r}\,\dot{\phi} + 2\,(r^2\sin\theta\cos\theta)\,\dot{\phi}\,\dot{\theta} \\\frac{\mathrm{d}}{\mathrm{d}\tau} \frac{\partial L}{\partial\dot{\theta}} &= r^2\,\ddot{\theta} + 2\,r\,\dot{r}\,\dot{\theta} \end{align}$$ And so $$\begin{align} \frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{r}}\right) - \frac{\partial L}{\partial r} &= \underbrace{1}_{g_{rr}}\,\ddot{r} + \underbrace{(-r\sin^2\theta)}_{\Gamma_{r\phi\phi}}\,\dot{\phi}^2 + \underbrace{(-r)}_{\Gamma_{r\theta\theta}}\,\dot{\theta}^2 \\\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{\phi}}\right) - \frac{\partial L}{\partial\phi} &= \underbrace{(r^2\sin^2\theta)}_{g_{\phi\phi}}\,\ddot{\phi} + 2\,\underbrace{(r\sin^2\theta)}_{\Gamma_{\phi r\phi} = \Gamma_{\phi\phi r}}\,\dot{r}\,\dot{\phi} + 2\,\underbrace{(r^2\sin\theta\cos\theta)}_{\Gamma_{\phi\phi\theta} = \Gamma_{\phi\theta\phi}}\,\dot{\phi}\,\dot{\theta} \\\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{\theta}}\right) - \frac{\partial L}{\partial\theta} &= \underbrace{r^2}_{g_{\theta\theta}}\,\ddot{\theta} + 2\,\underbrace{r}_{\Gamma_{\theta r\theta} = \Gamma_{\theta\theta r}}\,\dot{r}\,\dot{\theta} + \underbrace{(-r^2\sin\theta\cos\theta)}_{\Gamma_{\theta\phi\phi}} \,\dot{\phi}^2 \end{align}$$ All other Christoffel symbols of the first kind are zero. If we had a non-diagonal metric, some right-hand side expressions would have several second derivatives, each accompanied by a corresponding metric coefficient.

To obtain the Christoffel symbols of the second kind, find linear combinations of the above right-hand side expressions that leave only one second derivative, with coefficient $1$. Here this is easy because the metric is already in diagonal form. Therefore $$\begin{align} g^{rh}\left(\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{q}^h}\right) - \frac{\partial L}{\partial q^h}\right) &= \ddot{r} + \underbrace{(-r\sin^2\theta)}_{\Gamma^r_{\ \phi\phi}}\,\dot{\phi}^2 + \underbrace{(-r)}_{\Gamma^r_{\ \theta\theta}}\,\dot{\theta}^2 \\g^{\phi h}\left(\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{q}^h}\right) - \frac{\partial L}{\partial q^h}\right) &= \ddot{\phi} + 2\,\underbrace{\left(\frac{1}{r}\right)}_{\Gamma^\phi_{\ r\phi} = \Gamma^\phi_{\ \phi r}}\,\dot{r}\,\dot{\phi} + 2\,\underbrace{(\cot\theta)}_{\Gamma^\phi_{\ \phi\theta} = \Gamma^\phi_{\ \theta\phi}}\,\dot{\phi}\,\dot{\theta} \\g^{\theta h}\left(\frac{\mathrm{d}}{\mathrm{d}\tau} \left(\frac{\partial L}{\partial\dot{q}^h}\right) - \frac{\partial L}{\partial q^h}\right) &= \ddot{\theta} + 2\,\underbrace{\left(\frac{1}{r}\right)}_{\Gamma^\theta_{\ r\theta} = \Gamma^\theta_{\ \theta r}}\,\dot{r}\,\dot{\theta} + \underbrace{(-\sin\theta\cos\theta)}_{\Gamma^\theta_{\ \phi\phi}} \, \dot{\phi}^2 \end{align}$$ All other Christoffel symbols of the second kind are zero.

It's worth noting that if the metric is diagonal the calculations become much simpler. Explicitly,

\begin{align} &\Gamma^{\lambda}_{\;\,\mu \nu} = 0 \\ &\Gamma^{\lambda}_{\;\,\mu \mu} = -\frac{1}{2}(g_{\lambda \lambda})^{-1} \partial_{\lambda} g_{\mu \mu} \\ &\Gamma^{\lambda}_{\;\,\mu \lambda} = \partial_{\mu} (\ln \sqrt{|g_{\lambda \lambda}|}) \\ &\Gamma^{\lambda}_{\;\,\lambda \lambda} = \partial_{\lambda} (\ln \sqrt{|g_{\lambda \lambda}|}),\end{align}

where $\mu \neq \lambda \neq \nu \neq \mu$. These equations can easily be proven using the formula given in the original question.

You can read off the Christoffel symbols from the Euler–Lagrange equations for the Lagrangian naturally defined by the metric, $L=\frac12 g_{ij} v^i v^j$. Still messy, but at least there's not much you have to memorize in order to do the computation.

I can't recall a good reference right off the bat, but Google found this, for example: http://count.ucsc.edu/~rmont/classes/121A/polarChristoffelE_Lag.pdf