A system of $m$ (≥ 2) linear convection-diffusion two-point boundary value problems is examined, where the diffusion term in each equation is multiplied by a small parameter $ε$ and the equations are coupled through their convective and reactive terms via matrices $B$ and $A$ respectively. This system is in general singularly perturbed. Unlike the case of a single equation, it does not satisfy a conventional maximum principle. Certain hypotheses are placed on the coupling matrices $B$ and $A$ that ensure existence and uniqueness of a solution to the system and also permit boundary layers in the components of this solution at only one endpoint of the domain; these hypotheses can be regarded as a strong form of diagonal dominance of $B$. This solution is decomposed into a sum of regular and layer components. Bounds are established on these components and their derivatives to show explicitly their dependence on the small parameter $ε$. Finally, numerical methods consisting of upwinding on piecewise-uniform Shishkin meshes are proved to yield numerical solutions that are essentially first-order convergent, uniformly in $ε$, to the true solution in the discrete maximum norm. Numerical results on Shishkin meshes are presented to support these theoretical bounds.