Suppose $g$ is a Riemannian metric on a manifold, and $X$ is a smooth vector field. Show that $$(\mathcal{L}_Xg)_{ij} = \nabla_iX_j + \nabla_jX_i$$
I've been able to obtain that $$(\mathcal{L}_Xg)(\frac{\partial}{\partial x^i},\frac{\partial}{\partial x^j} ) = g(\nabla_{\frac{\partial}{\partial x^i}}X, \frac{\partial}{\partial x^j}) + g(\frac{\partial}{\partial x^i}, \nabla_{\frac{\partial}{\partial x^j}}X)$$
but when I expand this in local coordinates I get a bunch of pesky $g_{ij}$ terms that I can't make disappear. It also seems somewhat odd to me that the Lie derivative of a metric yields an answer that is a sum of derivatives of component functions of $X$. This makes it seem like the Lie Derivative of a metric is not dependent on the metric. Maybe I have misunderstood the author's notation, but I would like to see an explicit computation from my step to the final answer, as most sources I've found seem to gloss over it.
Your result is correct: one first computes $\mathcal L _X \big( g (Y,Z) \big)$, in this one replaces $X \cdot g(Y,Z)$ with the formula obtained from the fact that $\nabla$ is a metric connection, then using that $T(X,Y) = 0$ one gets $(\mathcal L _X) (Y, Z) = g(\nabla _Y X, Z) + g(Y, \nabla _Z X)$.
The author's formula is true with one condition, though: that instead of working with the local frame $\{ \partial _i \} _{i = 1, \dots, n}$ one should work with an orthonormal frame $\{ e _i \} _{i = 1, \dots, n}$. This follows easily from writing that $0 = \nabla _{e_i} g(e_j, e_k)$; assume that $\nabla _{e_i} e_a = A _{i a} ^b e_b$ (note that these are not necessarily Christoffel's coefficients) and do the computations having in mind that $g_{ij} = \delta _{ij}$ (Kronecker's delta). You will get that $A_{ji} ^k + A_{ki} ^j = 0$. Plug this in the result that you yourself have obtained and you will quickly obtain the author's formula.
To conclude, I suspect that the author must have made the assumption somewhere that he works with an orthonormal frame.