AbstractIn this article, we consider a system of convection-diffusion equations with delay terms. When a parameter multiplying the second order derivatives in the equations is small, boundary layers as well as interior layers appear in their solutions. A numerical method based on finite element and Shishkin mesh is presented. We derive an error estimate of orderO(N−1log2N) in the energy norm with respect to the parameter using the streamline-diffusion finite element method. Numerical experiments are also presented to support our theoretical results.