In various industrial combustion devices, such as liquid rocket engines at ignition or Diesel engines during the compression stage, the operating point varies over a wide range of pressures. These pressure variations can lead to a change of thermodynamic regime when the critical pressure is exceeded, switching from two-phase injection to transcritical injection. Such change modifies the topology of the flow and the mixing, thereby impacting the flame dynamics. This motivates the development of a unified methodology able to address both subcritical and supercritical flows within the same solver. To achieve this, the present work provides an extension of the supercritical real gas Taylor-Galerkin solver AVBP-RG to subcritical two-phase flows, based on diffuse interface models. In particular, the required developments for the integration of a multifluid model into the finite-element framework of this solver are detailed. Then, the ability of the solver to address a subcritical configuration is tested by simulating two subcritical-pressure operating points (G1 at 4.7 MPa and A10 at 1 MPa) of the MASCOTTE test bench operated by ONERA. This allows to confront the model with experimental data, showing good agreement.