  We consider a system of two singularly perturbed Boundary Value Problems
(BVPs) of convection-diffusion type with discontinuous source terms and a small
positive parameter multiplying the highest derivatives. Then their solutions
exhibit boundary layers as well as weak interior layers. A numerical method
based on finite element method (Shishkin and Bakhvalov-Shishkin meshes) is
presented. We derive an error estimate of order $O(N^{-1}\ln^{3/2}{N})$ in the
energy norm with respect to the perturbation parameter. Numerical experiments
are also presented to support our theoretical results.
