In this article, we extended operational matrices using orthonormal Boubaker polynomials of Riemann-Liouville fractional integration and Caputo derivative to find numerical solution of multi-term fractional-order differential equations (FDE). The proposed method is utilized to convert FDE into a system of algebraic equations. The convergence of the method is proved. Examples are given to explain the simplicity, computational time and accuracy of the method.