We consider a certain class of coupled systems of fractional differential equations involving ψ-Caputo fractional derivatives. A numerical approach is provided for solving this class of systems. The method is based on operational matrix of fractional integration of an arbitrary ψ-polynomial basis. A theoretical study related to the numerical scheme and the convergence of the method is presented. Next, several numerical examples are given using different types of polynomials aiming to confirm the efficiency of our approach.