Using the product rule for differentiation
$$ \nabla_{X} R(Y,Z)W + \nabla_{Y} R(Z,X)W + \nabla_{Z} R(X,Y)W
\\= \nabla_X(R(Y,Z)W) + \nabla_Y(R(Z,X)W) + \nabla_Z(R(X,Y)W)
\\- R(\nabla_X Y,Z)W - R(Y,\nabla_X Z)W
\\- R(\nabla_Y Z,X)W - R(Z,\nabla_Y X)W
\\- R(\nabla_Z X,Y)W - R(X,\nabla_Z Y)W
\\- R(Y,Z)\nabla_X W - R(Z,X)\nabla_Y W - R(X,Y)\nabla_Z W
$$
Also
$$R(\nabla_X Y,Z) + R(Y,\nabla_X Z) + R(\nabla_Y Z,X) + R(Z,\nabla_Y X) + R(\nabla_Z X,Y) + R(X,\nabla_Z Y)
\\= R(\nabla_X Y - \nabla_Y X,Z) + R(\nabla_Y Z - \nabla_Z Y,X) + R(\nabla_Z X - \nabla_X Z,Y)
\\= R(\tau(X,Y) + [X,Y],Z) + R(\tau(Y,Z) + [Y,Z],X) + R(\tau(Z,X) + [Z,X],Y)
$$
And
$$ \nabla_X(R(Y,Z)W) + \nabla_Y(R(Z,X)W) + \nabla_Z(R(X,Y)W)
\\= \nabla_X \nabla_Y \nabla_Z W - \nabla_X \nabla_Z \nabla_Y W
+ \nabla_Y \nabla_Z \nabla_X W - \nabla_Y \nabla_X \nabla_Z W
+ \nabla_Z \nabla_X \nabla_Y W - \nabla_Z \nabla_Y \nabla_X W
\\ - \nabla_X \nabla_{[Y,Z]} W - \nabla_Y \nabla_{[Z,X]} W - \nabla_Z \nabla_{[X,Y]} W
\\= \nabla_Y \nabla_Z \nabla_X W - \nabla_Z \nabla_Y \nabla_X W
+ \nabla_Z \nabla_X \nabla_Y W - \nabla_X \nabla_Z \nabla_Y W
+ \nabla_X \nabla_Y \nabla_Z W- \nabla_Y \nabla_X \nabla_Z W
\\ - \nabla_{[Y,Z]} \nabla_X W - R(X,[Y,Z]) - \nabla_{[X,[Y,Z]]}W - \nabla_{[Z,X]} \nabla_Y W - R(Y,[Z,X]) - \nabla_{[Y,[Z,X]]}W - \nabla_{[X,Y]} \nabla_Z W - R(Z,[X,Y]) - \nabla_{[Z,[X,Y]]}W
\\= R(Y,Z)\nabla_X W + R(Z,X)\nabla_Y W + R(X,Y)\nabla_Z W
\\ - R(X,[Y,Z])W - R(Y,[Z,X])W - R(Z,[X,Y])W
$$
noting that $\nabla_{[X,[Y,Z]]}W + \nabla_{[Y,[Z,X]]}W + \nabla_{[Z,[X,Y]]}W = \nabla_{[X,[Y,Z]] + [Y,[Z,X]] + [Z,[X,Y]]}W = 0$ by the Jacobi identity.
Now combine to get
$$ \nabla_{X} R(Y,Z) + \nabla_{Y} R(Z,X) + \nabla_{Z} R(X,Y)
= R(X,\tau(Y,Z)) + R(Y,\tau(Z,X)) + R(Z,\tau(X,Y))$$
Note the argument is simpler if one has that all the Lie brackets are zero, which is a much weaker assumption than normal coordinates, etc. I would say that the core of the proof is this identity:
$$ \nabla_X \nabla_Y \nabla_Z W - \nabla_X \nabla_Z \nabla_Y W
+ \nabla_Y \nabla_Z \nabla_X W - \nabla_Y \nabla_X \nabla_Z W
+ \nabla_Z \nabla_X \nabla_Y W - \nabla_Z \nabla_Y \nabla_X W
\\= \nabla_Y \nabla_Z \nabla_X W - \nabla_Z \nabla_Y \nabla_X W
+ \nabla_Z \nabla_X \nabla_Y W - \nabla_X \nabla_Z \nabla_Y W
+ \nabla_X \nabla_Y \nabla_Z W- \nabla_Y \nabla_X \nabla_Z W
$$
Finally let me comment that the problem with using normal coordinates is that you implicitly use that the connection comes from a Riemannian metric. The proofs given here are simple manipulations of the definitions.
The problem is, $e_ie_ie_j u$ and $e_je_ie_i u$ do not cancel. In fact, at $p$ we have
$e_ie_ie_j u-e_je_ie_i u=e_ie_ie_j u-e_ie_je_i u=e_i([e_i, e_j]u)$; now although $[e_i, e_j]=0$ at $p$, its derivative may not vanish; so $ e_i([e_i, e_j]u)$ may not be zero.
Now we compute
$(e_ie_ie_j u-e_je_ie_i u)(e_ju)=\Big(e_i([e_i, e_j]u)\Big)e_ju$; this is
$$e_i\langle [e_i, e_j], \nabla u\rangle e_ju
=e_i\langle \nabla_{e_i}e_j-\nabla_{e_j}e_i, \nabla u\rangle \cdot e_ju
=\langle \nabla_{e_i}\nabla_{e_i}e_j-\nabla_{e_i}\nabla_{e_j}e_i, \nabla u\rangle \cdot e_j u;
$$
the second term $-\langle \nabla_{e_i}\nabla_{e_j}e_i, \nabla u\rangle \cdot e_j u$ cancels the last remaining term you wrote; as for the first term, an easy way to get rid of it is, assume the frame is constructed like this: pick an orthonormal frame $e_i$ at $p$, then along any geodesic $\gamma$ from $p$ we parallel translate this frame along $\gamma$; this gives
a frame field locally defined near $p$. In particular, $\nabla_{e_i}e_j=0$ along the geodesic $\gamma_i$ at $p$ with initial tangent vector $e_i$, therefore $\nabla_{e_i}\nabla_{e_i} e_j=0$ along this geodesic $\gamma_i$, therefore at $p$ we have
$\nabla_{e_i}\nabla_{e_i} e_j=0$.
Best Answer
This is indeed a rather cumbersome calculation. Here it goes:
$$\nabla R(e_i,e_j,e_k,e_l,e_h)(p)+\nabla R(e_i,e_j,e_h,e_k,e_l)(p)+\nabla R(e_i,e_j,e_l,e_h,e_k)(p)=$$
$$\langle \nabla_{e_h} \nabla_{e_l} \nabla_{e_k}e_i - \nabla_{e_h} \nabla_{e_k} \nabla_{e_l}e_i +\nabla_{e_h} \nabla_{[e_k,e_l]}e_i$$
$$+\nabla_{e_l} \nabla_{e_k} \nabla_{e_h}e_i - \nabla_{e_l} \nabla_{e_h} \nabla_{e_k}e_i +\nabla_{e_l} \nabla_{[e_h,e_k]}e_i$$
$$+\nabla_{e_k} \nabla_{e_h} \nabla_{e_l}e_i - \nabla_{e_k} \nabla_{e_l} \nabla_{e_h}e_i +\nabla_{e_k} \nabla_{[e_l,e_h]}e_i, e_j \rangle (p) =$$
$$ \langle R(e_l,e_h)\nabla_{e_k}e_i-\nabla_{[e_l,e_h]}\nabla_{e_k}e_i+R(e_k,e_l) \nabla_{e_h}e_i -\nabla_{[e_k,e_l]}\nabla_{e_h}e_i+R(e_k,e_h)\nabla_{e_l}e_i-\nabla_{[e_k,e_h]}\nabla_{e_l}e_i+\nabla_{e_h} \nabla_{[e_k,e_l]}e_i+\nabla_{e_l} \nabla_{[e_h,e_k]}e_i+\nabla_{e_k} \nabla_{[e_l,e_h]}e_i,e_j\rangle (p)$$
$$= \langle R([e_l,e_h],e_k)e_i-\nabla_{[[e_l.e_h],e_k]}e_i+R([e_k,e_l],e_h)e_i-\nabla_{[[e_k.e_l],e_h]}e_i+R([e_k,e_h],e_l)e_i-\nabla_{[[e_k.e_h],e_l]}e_i,e_j\rangle(p) =$$
$$\langle R(\nabla_{e_l}e_h-\nabla_{e_h}e_l,e_k)e_i + R(\nabla_{e_k}e_l-\nabla_{e_l}e_k,e_h)e_i+R(\nabla_{e_k}e_h-\nabla_{e_h}e_k,e_l)e_i- \nabla_{[[e_l.e_h],e_k]+[[e_k.e_l],e_h]+[[e_k.e_h],e_l]}e_i,e_j\rangle(p) = 0$$
In this last line the Jacoby Identity and the symmetry of the connection is used. Observe that every time that elements of the form $\nabla_{e_l}e_h$ are dropped, what is being used is that of course $\nabla_{e_l}e_h(p)=0$, and the tensorial property of the Riemannian tensor (observe also the necessity that everything must be evaluated on p). Since p is arbitrary, and because of the linearity of everything involved, the identity is proven.