Abstract In this paper, we propose a new stabilized linear finite element method for solving reaction–convection–diffusion equations with arbitrary magnitudes of reaction and diffusion. The key feature of the new method is that the test function in the stabilization term is taken in the adjoint-operator-like form − e Δ v − ( a ⋅ ∇ v ) / γ + σ v , where the parameter γ is appropriately designed to adjust the convection strength to achieve high accuracy and stability. We derive the stability estimates for the finite element solutions and establish the explicit dependence of L 2 and H 1 error bounds on the diffusivity, modulus of the convection field, reaction coefficient and the mesh size. The analysis shows that the proposed method is suitable for a wide range of mesh Peclet numbers and mesh Damkohler numbers. More specifically, if the diffusivity e is sufficiently small with e ‖ a ‖ h and the reaction coefficient σ is large enough such that ‖ a ‖ σ h , then the method exhibits optimal convergence rates in both L 2 and H 1 norms. However, for a small reaction coefficient satisfying ‖ a ‖ ≥ σ h , the method behaves like the well-known streamline upwind/Petrov–Galerkin formulation of Brooks and Hughes. Several numerical examples exhibiting boundary or interior layers are given to demonstrate the high performance of the proposed method. Moreover, we apply the developed method to time-dependent reaction–convection–diffusion problems and simulation results show the efficiency of the approach.