We consider a system of singularly perturbed equations of reaction–diffusion type coupled through the reaction term. A high order Schwarz domain decomposition method is developed to solve the system numerically. The method splits the original domain into three overlapping subdomains. On two boundary layer subdomains we use a compact fourth order difference scheme on a uniform mesh while on the interior subdomain we use a hybrid scheme on a uniform mesh. We prove that the method is almost fourth order -uniformly convergent. Furthermore, we prove that when is small, one iteration is sufficient to get almost fourth order -uniform convergence. Numerical experiments are performed to support the theoretical results.