I am trying to plot a series molar volume versus pressure for Peng-robinson eos and volume translation equation, I am stuck here, please help me

Dears,

Thank you ahead for reading my questions and very appreciate your effort whether I can get solution from you or not.

I have tried a few nights and now I think I need another brilliant brain to help me out.

I am stuck in the part of %%%% v = Z*R*T(i)/P and can not run throught it.

My idea for this one is to solve Peng-robinson equation (with volume translation equation) for molar volume under a seriers of given pressure at isothermal process(different tempertaure isothermal), and finally plot the seriers of molar volume vs. pressure/reduced pressure at different tempertaure in a diagram.

Many thanks again.

Regards,

Please see my script as followis:

% CH4 is tested for supercritical phase  
clc
Tc = 190.546; % in K
Pc = 45.992; % in bar
omega = 0.011;
R = 0.083144; % l.bar/K/mol
a = 0.45724*R^2*Tc^2/Pc^2;
b = 0.07780*R*Tc/Pc;
kappa = 0.37464+1.54226*omega-0.26992*omega^2;
c1 = 0.01313; %component dependent parameter
% calculte conditions



for i = 1:20:100
    P=[46:46:460];
   
    T(i)= 190.546+i;
    Tr(i) = T(i)./Tc;
    Pr = P/Pc;
    alpha = (1+ kappa*(1-sqrt(Tr(i))))^2;
  
% calculate demensionless parameters
   A = alpha*a*P/(R*T(i))^2;
   B = b*P/R/T(i);

   % determine cubic coefficients and solve cubi
   m3 = 1;
   m2 = -(1-B);
   m1 = A-2*B-3*B.^2;
   m0 = -(A.*B-B.^2-B.^3);

% the following function finds the real roots.
% Zvals holds the real imaginary results.

   Zvals = roots([m3 m2 m1 m0]);

% determine indices for roots are real and greater than B
   index = find(imag(Zvals)==0);

% collect the real values of Z
   if length(index)>1
       %execute this if more than one root
       % collect the real values
      Zreal=real(Zvals(index));
      % accecpt the largeste and smallest real values because the center is
      % always unstable when three exist
      Z = [max(Zreal),min(Zreal)];
   else
       %execute this if there is just one real root
       Z = real(Zvals(index));
   end
       v = Z*R*T(i)/P;

d = v.^2/R/Tc*(R*T(i)./(v - b).^2)+ 2*a.*(v + b)./(v.^2 + 2*v.*b - b^2).^2;
delta = 0.092618044-0.105897139;
c = R*Tc/Pc*(c1-(0.004+c1)*exp(-2*d));
vvt= v + c -delta*(0.35./(0.35+d));


 h0=plot(v,P,'-g');
 hold on
 h=plot(vvt,P,'--r');


axis([0 1 0 500])
xlabel('Volume in l/mol')
ylabel('pressue in bar')
title ('Isotherms for Methane')
j=j+1;
end