Density calculation with cubic equation of state
First off, in the two-phase region you have three roots not two.

Smallest volume = liquid

Largest volume = gas

Middle volume = an unstable state you can't really reach

To determine whether a given temperature and pressure  gives the liquid or the gas you must compare the Gibbs free energy, $G$, between the liquid and gas roots as determined from solving the cubic equation.

The change in Gibbs free energy from liquid to gas is given by

$\Delta G = \int_{liquid}^{gas}{V dP}$

where $P$ and $V$ are pressure and volume along a path of constant temperature.  Note that the need for constant temperature; you can't just hold pressure constant and change the volume because that is not constant temperature.

The proper way to do it is to render

$dP=(dP/dV) dV$

where the derivative $dP/dV$ is obtained as a function of temperature and volume by differentiating the equation of state at constant temperature.  Then you have

$\Delta G = \int_{liquid}^{gas}{V(dP/dV) dV}$

with $dP/dV$ known as a function of volume at constant temperature.  You carry out the integral, which requires knowing how to integrate rational functions, and check whether $\Delta G$ is positive, negative, or zero:

$\Delta G&gt;0$ means liquid has lower free energy, so liquid is more stable.  You plot the point for the liquid phase root.

$\Delta G&lt;0$ means gas has lower free energy, so gas is more stable.  You plot the point for the gas phase root.

$\Delta G=0$ means both phases have the same free energy, you hit the equilibrium between liquid and gas. You plot both points and connect them with a line segment corresponding to both phases in equilibrium.  Jackpot!