Commutation of covariant derivative of functions Le$f$ be a smooth function on a Riemannian manifold $M$. My questions are:
a) If $\nabla_i f$ is a function, why is not true that $\nabla_j\nabla_k\nabla_if=\nabla_k\nabla_j\nabla_if$? This question arose when I wrote $\nabla_j\nabla_k\nabla_if=\nabla_j\nabla_k(\nabla_if)$, (is this true or not?). 
If we can see $\nabla_if$ as a function, then shouldn't be true that $\nabla_j\nabla_k(\nabla_if)=(\nabla^2f_i)(\partial_k,\partial_j)$, i. e., the Hessian of $\nabla_if$, which would give the symmetry on $k$ and $j$ on the formula $\nabla_j\nabla_k\nabla_if$? Where am I going wrong?
I know that there is something wrong on question a) because of the following formula:
$$
\nabla_j\nabla_i\nabla_jf-\nabla_i\nabla_j\nabla_jf=R_{jikj}\nabla_kf.
$$
b) How to prove the previous formula?
I appreciate any help.
 A: The reason why the covariant derivatives do not commute is precisely that they are not partial derivatives. If they were partial derivatives they would commute, but they are not.
For a function the covariant derivative is a partial derivative so 
$\nabla_i f = \partial_i f$
but what you obtain is now a vector field, and the covariant derivative, when it acts on a vector field has an extra term: the Christoffel symbol:
$\nabla_j (\partial_i f) = \partial_j\partial_i f + \Gamma_{ji}^m \partial_m f$
Due to the properties of the Christoffel symbols you can verify that in general the covariant derivatives will not commute with each other. For instance, above:
$\nabla_i \nabla_j f = \nabla_j \nabla_i f$
Only when there is no torsion, i.e. when $\Gamma^m_{ij} = \Gamma^m_{ji}$ (this needs not be true, but it is e.g. for the Levi-Civita connection).
In the particular example you gave you have to act with yet another covariant derivative $\nabla_k$, in this case when it acts on a tensor it will produce two Christoffel symbols. Putting all this together you may verify the equation you gave above for the curvature tensor $R_{ijkl}$
A: Your confusion seems to be notation related - in conventional abstract index notation, $\nabla_i (\nabla_j f)$ and $\nabla_i \nabla_j f$ are different things! The first means what you think it does - since $f$ is a function, $\nabla_j f = \partial_j f$ is also a function, and thus $\nabla_i(\nabla_j f) = \partial_i \partial_j f$. However, when the parentheses are omitted, all the covariant derivatives are taken before any indices are applied - that is, $$\nabla_i \nabla_j f := \nabla^2 f (\partial_i, \partial_j) = \nabla_i (\nabla_j f) - \nabla_{\nabla_i \partial_j} f.$$
This can be very confusing but it makes computations easier to write down.
Thus the resolution of your problem is that the expression
$$\nabla_j\nabla_k\nabla_if - \nabla_k\nabla_j\nabla_if$$is in fact commuting second covariant derivatives of a one-form $\nabla f$ - it is only at the very end that we take the $i$th component.
