A direct method for deriving an operational matrix of differentiation for Legendre polynomials is presented. The general case is also considered. The presented formulas are formally proved. Results are compared with the other known method. The new method is accurate for all matrix dimensions, while the other method may fail for higher dimensions.