In this article, we generalize the Legendre wavelets operational matrix of derivatives to fractional order derivatives in the Caputo sense. Legendre wavelets and their properties are employed for deriving Legendre wavelets operational matrix of fractional derivatives and a general procedure for forming this matrix is introduced. Then truncated Legendre wavelets expansions together with these matrices are used for numerical solution of Bagley–Torvik fractional order boundary value problems. Several examples are included to demonstrate accuracy and applicability of the proposed method. Key words: Shifted Legendre polynomials, Legendre wavelets, Caputo derivative, fractional order boundary value problems.