Philips in 1988 proved the following relationship:
If $f(x)$ is an infinitely differentiable function defined on the interval $[-1,+1]$ and its Legendre expansion is given by $f(x) = \sum_{n=0}^{\infty}a_{n}P_{n}(x)$, then the Legendre coefficients $a_{n}^{(q)}$ of the $q$-th derivative of $f(x)$ are given by $$a_{n}^{(q)} = \frac{(2n+1)}{2^{q-2}(q-1)!} \sum_{k=1}^{\infty} \frac{(k+q-2)!}{(k-1)!}\frac{(2n+2k+2q-3)!(n+k)!}{(2n+2k)!(n+q+k-2)!} a_{n+2k+q-2}$$.
Simplifying this for $q=1$, one finds
$$a_{n}^{(1)} = (2n+1) \sum_{k=1}^{\infty} a_{n+2k-1}$$
Please note that
You need a large number of coefficients (in theory an infinite amount) to calculate the coefficients of the derivative. This means that the order at which you cut-off the series will have an impact on the accuracy. If the original approximation is cut-off at order $M$, then clearly all $a_n$ for $n>M$ are identical to zero and the summation is finite.
Care must be taken in the summation if the $a_n$ have different signs. You'll quickly loose correct significant digits in IEEE double precision arithmetic.
The evaluation of $\sum_{n=0}^{N}a_{n}^{(q)}P_{n}(x)$ is best done using the Clenshaw-Smith algorithm (see Smith's paper). This is much more stable than converting to the monomial base and applying Horner's rule.
Reference
T.N. Philips, "On the Legendre Coefficients of a General-Order Derivative of an Infinitely Differential Function", IMA Journal of Numerical Analysis, 1988, Volume 8, p. 455--459