In this paper, the finite integration method and the operational matrix of fractional integration are implemented based on the shifted Chebyshev polynomial. They are utilized to devise two numerical procedures for solving the systems of fractional and classical integro-differential equations. The fractional derivatives are described in the Caputo sense. The devised procedure can be successfully applied to solve the stiff system of ODEs. To demonstrate the efficiency, accuracy and numerical convergence order of these procedures, several experimental examples are given. As a consequence, the numerical computations illustrate that our presented procedures achieve significant improvement in terms of accuracy with less computational cost.