In this paper, we have developed a new computational scheme for solving coupled systems of fractional order Volterra-type integro-differential equation (FVIDE). We construct new operational matrices, which serve as building blocks for converting the FVIDE into a Sylvester-type algebraic structure that is more easily solvable. The fractional derivatives and their inverses (integrals) are considered in the Caputo-Fabrizio sense. This computational scheme is based on Legendre polynomials. The assessment of the proposed method’s convergence involves utilizing appropriate error norms to measure the level of convergence. The algorithm is designed in such a way that it can be simulated using any computational package; we have used Matlab for simulating the proposed scheme. Graphical visualization and tabulation of data are performed to confirm convergence and to analyze errors.