Using a range of input to generate output, to plot. Using an Compressibility function.

Hi

I've a function written that calculates the Compressebility and Specific volume. I want to know how these output changes with changing Pressure(Which is an input for the function).

I want to work with a range of pressure indicated below, so i can eventually plot them.

P=0.9E06-P=9.5E06

My input:

Pc=[4.599E06 4.872E06 4.248E06 3.796E06 3.37E06 3.025E06 2.49E06 1.817E06 1.401E06]; 
Tc=[190.6 305.3 369.8 425.1 469.7 507.6 568.7 658.1 722];
T=320;
w=[0.012 0.1 0.152 0.2 0.252 0.301 0.346 0.577 0.721];
P=E06;
kappa=0;
eta=0;
x=[0.7 0.08 0.07 0.03 0.01 0.02 0.04 0.02 0.03];

My main Function:

function [Z,V,ZL,ZV,VL,VV] = ZMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x)
%The funtion ZMIX_PRVDW calculates the compressibility factor z and the specific
%molar volume V for a mixture at given pressure P, temperature T and concentration
%using PR-CEoS and VDW MIXING RULES. This function requires REALROOTS_CARDANO and
% MIXRULES_VDW funtions to run.
%All input/output data are expressed in SI.
%P[Pa], T[K], w[dimensionless], V[m^3/mol], Z[dimensionless]

%universal gas constant
R=8.314;

%PR attractive constant and covolume at critical point
b=0.07780*(R*Tc)./Pc;
k=0.37464+1.54226*w-0.26992*w.^2;
alpha=(1+k.*(1-sqrt(T./Tc))).^2;
a=0.45724*alpha.*(R^2.*(Tc.^2))./Pc;
c=length(x);

%calling MIXRULES_VDW function
aa=zeros(c,c);
bb=zeros(c,c); 
[am, bm, aa] = MIXRULES_VDW(a, b, kappa, eta, x);

%calculation of Am and Bm
Am=am*P/(R*T)^2;
Bm=bm*P/(R*T);

%calculation of A and B
a2=-1+Bm;
a1=Am-3*Bm.^2-2*Bm;
a0=-Am.*Bm+Bm.^2+Bm.^3;

%calculation of Z and V using REALROOTS_CARDANO function
Z=REALROOTS_CARDANO(a2,a1,a0);
V=Z*R*T/P;

%selection of Z and V values for liquid and vapor
if min(V)<b
    ZL=max(Z);
    ZV=max(Z);
    VL=max(V);
    VV=max(V);
else
    ZL=min(Z);
    ZV=max(Z);
    VL=min(V);
    VV=max(V);
    
end

end

My Subprograms:

function [am, bm, aa] = MIXRULES_VDW(a, b, kappa, eta, x)
%The funtion MIXRULES_VDW calculates the am and bm of a mixture using  
%VDW MIXING RULES. 

%All input/output data are expressed in SI or are dimensionless.

%giving a firt value for c
c=length(x);

%initialization the values of am and bm
am=0;
bm=0;

%calculation of am, bm and aa
aa=zeros(c,c); 
bb= zeros(c,c); 
aa = sqrt(a .* a.') .* (1-kappa);  
bb = (b + b.')/2 .* (1-eta);       
am = sum(x .* x.' .* aa, 'all');   
bm = sum(x .* x.' .* bb, 'all');

Real root subprograms used to find roots:

function [RR]=REALROOTS_CARDANO(a2,a1,a0)
%This function calculates the real roots of a polynomial of degree 3 using 
%an algorithm based on the Cardano's formula. 
%The polynomial must be written as: x^3 + a2*x^2 + a1*x + a0. 
%Input data: a2, a1, a0.
%Output data: array containing the real roots. 
%The output array can be 1x1 or 1x3. 

Q=(a2.^2-3*a1)/9;
R=(2*a2.^3-9*a1*a2+27*a0)/54;

    if Q==0 && R==0 %Case of single real root with multiplicity 3
        RR=-a2/3;
    else
        if Q^3-R^2>=0 
            theta=acos(R/(Q^(3/2)));
            RR(1)=-2*sqrt(Q)*cos(theta/3)-a2/3;
            RR(2)=-2*sqrt(Q)*cos((theta+2*pi)/3)-a2/3;
            RR(3)=-2*sqrt(Q)*cos((theta+4*pi)/3)-a2/3;
        else
            RR=-sign(R)*((sqrt(R^2-Q^3)+abs(R))^(1/3)+Q/(sqrt(R^2-Q^3)+abs(R))^(1/3))-a2/3;
        end
    end
end