Clear["Global`*"]

toSI[x_] := N[QuantityMagnitude[UnitConvert[Quantity[1, x], "SIBase"]]]

p[V_, T_] := ((n k T)/(V - n b))- ((a (n^2))/(V^2))

G[V_, T_] := - n k T Log[V - n b] + (((n k T) ( n b))/(V - n b)) - ((2 a (n^2))/V)

Vc = 3 n b;

n = toSI["AvogadroConstant"]

k = toSI["BoltzmannConstant"]

b = 6. * (10^(-29)); a = 4. * (10^(-49));

Vc = 3 n b

pc = (a/(27 (b^2)))

Tc = ((8 a)/(k 27 b))

Plot[Evaluate[Table[p[V, Ts], {Ts, Tc-40, Tc +40, 10}]], {V, Vc/2, 2 Vc}, PlotRange -> All, AxesLabel->Larger]

Plot[Evaluate[Table[G[V, Ts], {Ts, Tc-40, Tc +40, 10}]], {V, Vc/2, 3 Vc}, PlotRange -> All, AxesLabel->Larger]

data = Table[
Select[Table[
{p[V, Tc+dT], G[V, Tc+dT]}, {V, Vc/4, 20 Vc, Vc/4}], Head[#[[2]]] == Real &],
{dT, -20, 8, 4}];

ListLinePlot[data, LabelStyle->Larger, AxesLabel -> {"Pressure (Pa)", "Gibbs Free Energy (J)"}, ImageSize -> Scaled[0.8]]