However, I had difficulty understanding that answer and would like to understand how to do it this way. That is to say, I'd really like to know what property or identity that I'm missing before I can use use the Bianchi identities to show that it is manifestly zero.
The other proof uses the first Bianchi identity. That's where the starting assumption $R^a{}_{bcd}\xi^d = \xi^a{}_{;bc}$ comes from. If you want to use the second Bianchi identity, it is
$$(\nabla_\xi R)(X,Y) + (\nabla_X R)(Y,\xi) + (\nabla_Y R)(\xi,X) = 0\text{,}$$
and therefore applying it and the Leibniz rule produces:
$$\begin{align}
\underbrace{\nabla_\xi[R(X,Y)]+\nabla_X[R(Y,\xi)]+\nabla_Y[R(\xi,X)]}_\mathrm{foo} = \underbrace{R(\mathcal{L}_\xi X,Y) + R(\mathcal{L}_XY,\xi) + R(\mathcal{L}_Y\xi,X)}_\mathrm{bar}\text{,}
\end{align}$$
where it was assumed that the torsion vanishes, so that $\mathcal{L}_AB = \nabla_AB-\nabla_BA$. Additionally,
$$\begin{eqnarray*}
(\mathcal{L}_\xi R)(X,Y) &=& \mathcal{L}_\xi[R(X,Y)] - R(\mathcal{L}_\xi X,Y) - R(X,\mathcal{L}_\xi Y)\\
&=&\mathcal{L}_\xi[R(X,Y)] - R(\mathcal{L}_\xi X,Y) - R(\mathcal{L}_Y\xi,X)\\
&=&\underbrace{\mathcal{L}_\xi[R(X,Y)] - [\mathrm{foo}] + R(\mathcal{L}_XY,\xi)}_\mathrm{qux}\text{.}
\end{eqnarray*}$$
So the objective is to show that the right-hand side, $\mathrm{qux}$, is identically zero whenever $\xi$ is a Killing vector field.
Let's write $S^a{}_b = [R(X,Y)]^a{}_b = R^a{}_{bcd}X^cY^d$, and just crank it out:
$$\begin{eqnarray*}
\mathcal{L}_\xi S^a{}_b &=& \nabla_\xi S^a{}_b - S^e{}_b\xi^a{}_{;e} + S^a{}_e\xi^e{}_{;b}\\
&=& \nabla_\xi S^a{}_b + X^cY^d(R^a{}_{ecd}\xi^e{}_{;b} - R^e{}_{bcd}\xi^a{}_{;e})\\
&=& \nabla_\xi S^a{}_b + X^cY^d(\nabla_c\nabla_d-\nabla_d\nabla_c)\xi^a{}_{;b}\text{,}
\end{eqnarray*}$$
where the last step is actually valid for arbitrary $Z^a{}_b$, not just $\xi^a{}_{;b}$. The first term of this cancels with the first term of $\mathrm{foo}$. So far we have not used the fact that $\xi$ is a Killing vector field. Let's do so now by considering the other two terms of $\mathrm{foo}$:
$$\nabla_X[R(Y,\xi)]^a{}_b - \nabla_Y[R(X,\xi)]^a{}_b = \nabla_X\nabla_Y\xi^a{}_{;b} - \nabla_Y\nabla_X\xi^a{}_{;b}\text{,}$$
where the starting identity $R^a{}_{bcd}\xi^d = \xi^a{}_{;bc}$ was used. The same identity also gives:
$$R(\mathcal{L}_XY,\xi)^a{}_{b} = \nabla_{[X,Y]}\xi^a{}_{;b}\text{.}$$
Therefore, we have shown that for any vector fields $X,Y$,
$$\begin{eqnarray*}
X^cY^d(\mathcal{L}_\xi R^a{}_{bcd}) &=& \left[X^cY^d(\nabla_c\nabla_d-\nabla_d\nabla_c)-(\nabla_X\nabla_Y-\nabla_Y\nabla_X) + \nabla_{[X,Y]}\right]\xi^a{}_{;b}\\
&=& 0\text{.}\end{eqnarray*}$$
(If you have trouble with the last step, check Christoph's answer to the other question and modify appropriately.) Thus $\mathcal{L}_\xi R^a{}_{bcd} = 0$, QED.
One defining property of Christoffel symbols of the second kind is
$d\mathbf{e}_i=\Gamma^k_{ij}\mathbf{e}_k dq^j$.
Accepting this as a definition for the object $\Gamma^k_{ij}$ one can show, looking at the second derivative of the line element, that $\Gamma$ is symmetrical in its lower indices $\Gamma^k_{ij}=\Gamma^k_{ji}$.
Now to the derivation of an expression for $\Gamma$: looking at the total derivative of the metric one can get to:
$dg_{ij}=d( \mathbf{e}_i \cdot \mathbf{e}_j )=(\Gamma^k_{jl}g_{ik}+\Gamma^k_{il}g_{jk})dq^l$.
But by definition the total derivativ of $g_{ij}$ is given by $dg_{ij}=\frac{\partial g_{ij}}{\partial q^l}dq^l$. By compare the coefficients we get to:
$\frac{\partial g_{ij}}{\partial q^l}=\Gamma^k_{jl}g_{ik}+\Gamma^k_{il}g_{jk}$.
EDIT: That is what was derived in the question by a different way. But now to isolate a single Christoffel symbol one needs to add this expression up with different indicies. The mistake in the derivation in the question was pointed out in the comments; it was a mistake concering the summation index $\lambda$.
Using that one can show using the symmetries of $\Gamma$ and $g$ that the following holds:
$\frac{\partial g_{ij}}{\partial q^l}+\frac{\partial g_{lj}}{\partial q^i}-\frac{\partial g_{il}}{\partial q^j}=2\Gamma^k_{li}g_{jk} $.
Now dividing by 2 and inverting with $g$ gets you to an expression for $\Gamma$:
$\Gamma^k_{li}=\frac{1}{2}g^{jm}(\frac{\partial g_{ij}}{\partial q^l}+\frac{\partial g_{lj}}{\partial q^i}-\frac{\partial g_{il}}{\partial q^j})$
This is one possible derivation where granted the step of summing up those 3 partial derivatives is not very intuitive.
I know one can get to an expression for the Christoffel symbols of the second kind by looking at the Lagrange equation of motion for a free particle on a curved surface. This basically get you the geodetic equation where $\Gamma$ shows up as well. Then the devining property would be the geodetic equation and one would need to do the above calculation to show that
$d\mathbf{e}_i=\Gamma^k_{ij}\mathbf{e}_k dq^j$ actually holds for $\Gamma$.
Best Answer
The covariant derivative for a general tensor of the form $T^{a_1\dots a_n}_{b_1 \dots b_n}$ is given by,
$$\nabla_c T^{a_1\dots a_n}_{b_1 \dots b_n} = \partial_c T^{a_1\dots a_n}_{b_1 \dots b_n} + \Gamma^{a_1}_{cd}T^{d\dots a_n}_{b_1 \dots b_n} + \dots - \Gamma^d_{c b_1}T^{a_1\dots a_n}_{d \dots b_n} - \dots$$
Taking the covariant derivative of a covariant field $V_a$, we find,
$$\nabla_b V_a = \partial_b V_a - \Gamma^c_{ba}V_c$$
Now, the object $\nabla_b V_a$ has two lower indices, so taking the covariant derivative again, we find,
$$\nabla_c (\nabla_b V_a) = \partial_c(\nabla_b V_a) - \Gamma^d_{cb} (\nabla_d V_a) - \Gamma^d_{ca}(\nabla_b V_d)$$
Inserting the original covariant derivative, we find explicitly,
$$\nabla_c (\nabla_b V_a) = \partial_c (\partial_b V_a -\Gamma^{e}_{ba}V_e) - \Gamma^d_{cb}(\partial_d V_a - \Gamma^e_{da}V_e) - \Gamma^d_{ca}(\partial_b V_d - \Gamma^e_{bd}V_e)$$