In this article, a fully discrete numerical method is studied to solve the system of singularly perturbed parabolic convection–diffusion problem. The full discretization incorporates the backward-Euler method and the streamline-diffusion finite element method (SDFEM) for the temporal and spatial discretization, respectively. In each equation, the highest ordered spatial derivatives are multiplied by distinct parameters, hence the overlapping boundary layers phenomena displayed in the solution. To capture the singularity of the solution, a layer-adapted mesh has been used for the spatial domain discretization. The stability of the numerical method has been addressed by considering an appropriate stabilization parameter. A complete error analysis has shown the first order of accuracy of the fully discrete method in the discrete -norm. Lastly, few numerical experiments have been addressed to justify the theoretical estimates.