A parameter-uniform numerical scheme for a system of weakly coupled singularly perturbed reaction-diffusion equations of arbitrary size with appropriate boundary conditions is investigated. More precisely, quadratic $B$-spline basis functions with an exponentially graded mesh are used to solve a $\ell\times\ell$ system whose solution exhibits parabolic (or exponential) boundary layers at both endpoints of the domain. A suitable mesh generating function is used to generate the exponentially graded mesh. The decomposition of the solution into regular and singular components is obtained to provide error estimates. A convergence analysis is addressed, which shows a uniform convergence of the second order. To validate the theoretical findings, two test problems are solved numerically.