The numerical simulation of full-disequilibrium multi-phase flows is intrinsically expensive due to the number of equations involved. The addition of stiff mechanical relaxation and phase transition increases the cost further. For these reasons, most work has been limited to numerical methods using simple equations of state, such as the Noble-Able Stiffened Gas. When simulating molecularly complex fluids near saturation or critical conditions, more accurate equations of state, such as the reference equation of state by Span and Wagner [1], are necessary to capture the correct thermodynamic behavior. Because these equations of state are naturally expressed in terms of density and temperature, their use in a flow solver, which generally provides only density and internal energy directly, can be extremely expensive. To accurately model fluid flows near saturation, we present a finite-volume second-order solver for the Baer-Nunziato model, paired with the thermodynamic library CoolProp [2], that minimizes the cost associated with the evaluation of the arbitrary equation of state thanks to a conservative temperature update of the hy- perbolic part of the phasic total-energy equations. We also discuss a modified version of the stiff mechanical relaxation solver of Chiocchetti and M¨uller [3] for arbitrary equations of state. Multiple tests considering different fluids show that the method successfully described two-phase flows evolving close to the saturation curve and experiencing phase transition, thanks to a robust strategy to handle the instabilities in the equation of state inside the saturation curve.