Exact equations

v, p = var('v p')
palette = [(215/255, 0/255, 132/255), (255/255, 1/255, 73/255), (255/255, 121/255, 1/255), (255/255, 210/255, 0/255)]
cool_palette = [(0/255, 150/255, 173/255), (0/255, 200/255, 146/255)]

The van der Waals equation

energy_densities = [k/27 for k in range(-45, 54)]
const_energy_curves = [plot(en/(v-1) - 1/v^2, (v, 1, 6), ymin = -1/3, ymax = 1/3, rgbcolor = cool_palette[0]) for en in energy_densities]
gas_labels = ["$v$", "$p$"]
show(sum(const_energy_curves) + point((3, 1/27), color = palette[2], size = 16), axes_labels = gas_labels)