The aim of this paper to investigate a weak Galerkin finite element method (WG-FEM) for solving a system of coupled singularly perturbed reaction-diffusion equations. Each equation in the system has perturbation parameter of different magnitude and thus, the solutions will exhibit two distinct but overlapping boundary layers near each boundary of the domain. The proposed method is applied to the coupled system on Shishkin mesh to solve the problem theoretically and numerically. Elimination of the interior unknowns efficiently from the discrete solution system reduces the degrees of freedom and, thus the number of unknown in the discrete solution is comparable with the standard finite element scheme. The stability and error analysis of the proposed method on the Shishkin mesh are presented. We show that the method convergences of order O(N −k ln k N) in the energy norm, uniformly with respect to the perturbation parameter. Moreover, the optimal convergence rate of O(N −(k+1)) in the L 2-norm and the superconvergence rate of O((N −2k ln 2k N) in the discrete L ∞-norm is observed numerically. Finally, some numerical experiments are carried out to verify numerically theory.