<p>This study introduces a new and efficient way to calculate operational matrices for Caputo and modified Atangana Baleanu fractional derivatives in the framework of generalized Morgan Voyce polynomials. These matrices serve as the foundation for devising a computational framework to address an abstract multi-order fractional integro-differential equation, encompassing well-established cases such as the Bagley–Torvik, Duffing, and Painlev’e equations, among others. These equations have significant applications in science and engineering, such as in fluid dynamics, modelling of viscoelastic materials, oscillatory systems with memory, nonlinear phenomena in physics and optical systems. Initially, the integer and fractional order operational matrices are introduced to transform the problem and its constraints into a system of algebraic equations with unknown coefficients, which are used to approximate the solution. Further, the collocation strategy is employed to solve the system of algebraic equations. Thus, these coefficients facilitate determining solutions for the proposed model. This methodology offers a straightforward implementation and exhibits enhanced computational efficiency relative to conventional techniques, requiring a reduced number of basis functions and demanding less computational effort. To evaluate the effectiveness of the proposed method, several examples are simulated, and the results are compared with both analytical solutions and previously reported findings in the literature.</p>