The present work is devoted to developing two numerical techniques based on fractional Bernstein polynomials, namely fractional Bernstein operational matrix method, to numerically solving a class of fractional integro-differential equations (FIDEs). The first scheme is introduced based on the idea of operational matrices generated using integration, whereas the second one is based on operational matrices of differentiation using the collocation technique. We apply the Riemann–Liouville and fractional derivative in Caputo’s sense on Bernstein polynomials, to obtain the approximate solutions of the proposed FIDEs. We also provide the residual correction procedure for both methods to estimate the absolute errors. Some results of the perturbation and stability analysis of the methods are theoretically and practically presented. We demonstrate the applicability and accuracy of the proposed schemes by a series of numerical examples. The numerical simulation exactly meets the exact solution and reaches almost zero absolute error whenever the exact solution is a polynomial. We compare the algorithms with some known analytic and semi-analytic methods. As a result, our algorithm based on the Bernstein series solution methods yield better results and show outstanding and optimal performance with high accuracy orders compared with those obtained from the optimal homotopy asymptotic method, standard and perturbed least squares method, CAS and Legendre wavelets method, and fractional Euler wavelet method.