<p>This paper analyses the analytical approximation solution for linear and nonlinear systems of fractional differential equations (FDEs) using shifted Legendre polynomials. We show how to solve nonlinear fractional problems by applying the operational matrix technique with the Caputo derivative. The Legendre operational matrix, derived using orthogonality properties, is implemented in dynamical fractional problems. By transforming the FDEs into a system of nonlinear algebraic equations, and applying Newton’s iterative scheme with an initial guess, we obtain polynomial coefficients that yield the approximate solution. We provide convergence analysis and an illustrative example to validate the proposed method. Furthermore, we analyse the nonlinear fractional SIS epidemic model and nonlinear coupled systems using the spectral collocation method. Numerical simulations and graphical results demonstrate the efficacy of the spectral collocation and Euler methods for solving real-world problems involving fractional-order differential equations.</p>