Generalise $\nabla_Y(fX)=Y(f)X+f\nabla_Y(X)$ First an example: suppose we are on a $2$-sphere, and we have two vector fields $X,Y$ on it. Then, suppose that we want to rotate all vectors of $X$ by $\frac{\pi}{2}$ with respect to the outer normal:
in spherical local coordinates ($\theta\in (0,\pi)$, $\varphi\in (0,2\pi)$), this rotation corresponds to the matrix $R(\theta)=\begin{bmatrix}0 & -\sin(\theta)\\\frac{1}{\sin(\theta)}& 0\end{bmatrix}$. We would like to compute $\nabla_Y(R(\theta)X)$.
My question: is there any property of the covariant derivative that allows me to compute easily $\nabla_Y(A(x)X)$, where $A$ is some coordinates dependant matrix acting in local coordinates on $Y$?

Something that generalise $\nabla_Y(fX)=Y(f)X+f\nabla_Y(X)$ to matrices.

Here the problem is that the matrix actually mix together the vector components $\frac{\partial}{\partial X_i}$.
 A: To answer your question about rotation: First, the rotation $R$ by $\pi/2$ along the normal vector $\vec n$ is given by 
$$ R (X) =\vec n \times X,$$
where $\times $ is the cross product in $\mathbb R^3$. For any other $Y$, 
\begin{align}
Y( R(X) ) &= \nabla_Y (\vec n \times X) \\
&= Y(\vec n) \times X + \vec n \times Y(X).
\end{align}
By $Y(R(X))$ and $Y(X)$ I mean the componentwise differentiation. Note that $\nabla _YX = (Y(X))^T$ by definition of $\nabla$. For the $2$-sphere, $x\mapsto \vec n(x)$ is indeed the immersion $\mathbb S^2 \subset \mathbb R^3$. Thus $Y(\vec n) = Y$. Also, $\vec n \times Y(X) = \vec n \times \nabla_Y X$, so 
$$Y(R(X)) = Y\times X + R(\nabla_YX)\Rightarrow \nabla_Y(R(X)) = R(\nabla_YX).$$
We can do something similar for any immersed surface $M$ in $\mathbb R^3$. We still have $R(X) = \vec n \times X$ and 
$$ Y(R(X)) = Y(\vec n) \times X + \vec n \times \nabla_YX. $$
This time $Y(\vec n) = A(Y)$, where $A$ is the shape operator, which maps tangent vectors to tangent vectors. In particular, $A(Y) \times X$ is orthogonal to your surface, thus we also have 
$$ \nabla _Y(R(X)) = R(\nabla_YX).$$
In general, for any $\theta$, the rotation $R(\theta)$ is given by 
$$R(\theta) X = \cos \theta X + \sin \theta R(X).$$
Thus 
\begin{align}
\nabla_Y ( R(\theta)(X)) &= \cos \theta \nabla _Y(X) + \sin\theta \nabla_Y( R(X)) \\
&=  \cos \theta \nabla _Y(X) + \sin\theta R(\nabla_YX) \\
&= R(\theta) (\nabla_YX). 
\end{align}
Thus for any immersed surface $M$ in $\mathbb R^3$, the rotation $R(\theta)$ commutes with covariant derivatives. 
A: In general, for an endomorphism field $A$ of $TM$, any covariant derivative satisfies the Leibniz Rule
$$\nabla_Y (AX) = (\nabla_Y A) X + A \nabla_Y X .$$
(But NB often this is how $\nabla_Y A$ is defined in the first place.)
In a coordinate frame with Christoffel symbols $\Gamma^d_{ef}$,
$$\nabla_c A^a{}_b = \partial_c A^a{}_b + \Gamma^a_{dc} A^d{}_b - \Gamma^d_{bc} A^a{}_d .$$
With respect to any frame, we can interpret $$\nabla_c A^a{}_b$$ as a matrix of $1$-forms indexed by $(a, b)$.
In our particular example of the round $2$-sphere, using the round metric $g$ to raise the second index of the usual volume form $\operatorname{vol}$ produces an endomorphism field $J$, and one can check that it coincides with the map $R$ that rotates each tangent space $\frac{\pi}{2}$ w.r.t. the outward-pointing normal.
But the Levi-Civita connection $\nabla$ satisfies $\nabla g = 0$ by definition, and since (once we've fixed an orientation) $\operatorname{vol}$ is defined invariantly using $g$, we also have $\nabla \operatorname{vol} = 0$ and thus $\nabla J = 0$. (If you prefer a more concrete proof, we can check this in local coordinates, too.) So, our Leibniz rule gives
$$\nabla_Y(JX) = J \nabla_Y X .$$
But from the wording of the original question I believe you already know how to evaluate the right-hand side.
Remark The same proof works just as well for any oriented $2$-manifold. (Reversing the choice of orientation replaces $J$ with its negative.)
