I don't know what you are looking exactly. I did something here, please let me know is it correct or you need further help
clc; clear all; close all;
T = 310:5:350;
Tc = 369.8;
Tr1 = T./Tc;
Vr = linspace(0.5,5,500);
figure
for ii = 1:length(Tr1)
Tr = Tr1(ii);
funcPr = @(Vr) 8/3 * Tr./(Vr - 1/3) - 3./(Vr.^2);
Pr = funcPr(Vr);
if Tr < 1
Pr_l = 1.0;
vdW_Pr_l = [1 -1/3*(1+8*Tr/Pr_l) 3/Pr_l -1/Pr_l];
v=sort(roots(vdW_Pr_l));
A1 = (v(2)-v(1))*Pr_l - integral(funcPr,v(1),v(2));
A2 = integral(funcPr,v(2),v(3)) - (v(3)-v(2))*Pr_l;
AreaDiff = abs (A1-A2);
while AreaDiff > 0.0001
vdW_Pr_l = [1 -1/3*(1+8*Tr/Pr_l) 3/Pr_l -1/Pr_l];
v=sort(roots(vdW_Pr_l));
funcPr = @(Vr) 8/3 * Tr./(Vr - 1/3) - 3./(Vr.^2);
A1 = (v(2)-v(1))*Pr_l - integral(funcPr,v(1),v(2));
A2 = integral(funcPr,v(2),v(3)) - (v(3)-v(2))*Pr_l;
AreaDiff = abs (A1-A2);
Pr_l = Pr_l - 0.00001;
end
end
plot(Vr,Pr); hold all
xlim([0.25 4])
ylim([-0.5 2])
xlabel('Vr')
ylabel('Pr')
plot([0.25 4],[Pr_l Pr_l],'k--')
plot(v(1),Pr_l,'bo',v(2),Pr_l,'bo',v(3),Pr_l,'bo')
end
