In the present article, an effective Laguerre collocation method is used to obtain the approximate solution of a system of coupled fractional order non-linear reaction-advection-diffusion equations with prescribed initial and boundary conditions. In the proposed scheme, Laguerre polynomials are used together with operational matrix and collocation method to obtain approximate solutions of the coupled systems, so that our proposed model is converted into a system of algebraic equations which can be solved employing the Newton method. The solution profiles of the coupled systems are presented graphically for different particular cases. The salient features of the present article are finding the stability analysis of the proposed method and also the demonstration of the lower variations of solute concentrations with respect to the column length in fractional order system compared to integer order system. To show the higher efficiency, reliability and accuracy of the proposed scheme, a comparison between the numerical results of Burgers’ coupled systems and its existing analytical result is reported. There are high compatibility and consistency between the approximate solution and its exact solution to a higher order of accuracy. The exhibition of error analysis for each case through tables and graphs confirms the super-linearly convergence rate of the proposed method.