T=200; %K
a=1.505; %L^2*bar/mol^2
b=0.03985; %L/mol
R=8.3144598*10^-2; %L*?bar/K*?mol

for P=(1:1:5) %range on x axis
    A=(a.*P)/((R*T)^2)
    B=(b.*P)/(R*T)
    r=roots([1,-(1+B),A,-A.*B]); %finding roots
    Z(P)=r(imag(r)==0); %keeping only real roots
    gamma=exp((Z-1)-(A./Z)-(log(Z-B))) %final equation
end
F=(1:1:5); %range once again

plot(F,gamma)