We show that the variational energy principle of the multi-region relaxed magnetohydrodynamic (MRxMHD) model can be used to predict finite-pressure linear tearing instabilities. In this model, the plasma volume is sliced into sub-volumes separated by 'ideal interfaces', and in each volume the magnetic field relaxes to a Taylor state, where the pressure gradient $\nabla p = 0$. The MRxMHD model is implemented in the Stepped-Pressure Equilibrium Code (SPEC) so that the equilibrium solution in each region is computed while preserving the force balance across the interfaces. As SPEC computes the Hessian matrix (a discretized stability matrix), the stability of an MRxMHD equilibrium can also be computed with SPEC. In this article, using SPEC, we investigate the effect of local pressure gradients and the $\nabla p = 0$ in the vicinity of the resonant surface of a tearing mode. For low-beta plasma, we have been able to illustrate a relationship between the resistive singular-layer theory (Coppi et al 1966 Nucl. Fusion6 101; Glasser et al 1975 Phys. Fluids18 875–88) and the MRxMHD model. Within the singular layer, the volume-averaged magnetic helicity and the flux-averaged toroidal flux are shown to be the invariants for the linear tearing modes in SPEC simulations. Our technique to compute MRxMHD stability is first tested numerically in a cylindrical tokamak and its application in toroidal geometry is demonstrated. We demonstrate an agreement between the stability boundary obtained with SPEC simulation and the resistive inner-layer theories.