In this paper, we generalize a coupled system of nonlinear reaction-advection-diffusion equations to a variable-order fractional one by using the Caputo-Fabrizio fractional derivative, which is a non-singular fractional derivative operator. In order to establish an appropriate method for this system, we introduce a new formulation of the discrete Legendre polynomials namely the orthonormal shifted discrete Legendre polynomials. The operational matrices of classical and fractional derivatives of these basis functions are extracted. The devised method uses these polynomials and their operational matrices together with the collocation technique to transform the system under consideration into a system of algebraic equations which is uncomplicated for solving. Two numerical examples are analyzed to examine the accuracy of the method.