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.　©　2020