Integro-differential equations are developed as models in enormous fields of engineering and science such as biological models, population growth, aerospace systems and industrial mathematics. In this work, we consider a general class of nonlinear fractional integro-differential equations with variable order derivative. We use the operational matrices based on the Bernstein polynomials to obtain numerical solution of this type of equations. By utilizing the operational matrices along with the Newton–Cotes collocation points, the problem under study is reduced to a system of nonlinear algebraic equations. An error estimate of the numerical solution is proved. Finally, some examples are included to show the accuracy and validity of the proposed method.