This paper is devoted to the approximation of the convection–diffusion–reaction equation using a mixed, first-order, formulation. We propose, and analyse, a stabilised finite element method that allows equal order interpolations for the primal and dual variables. This formulation, reminiscent of the Galerkin least-squares method, is proven stable and convergent. In addition, a numerical assessment of the numerical performance of different stabilised finite element methods for the mixed formulation is carried out, and the different methods are compared in terms of accuracy, stability, and sharpness of the layers for two different classical test problems.