One point calculation to a range of points.

clear all
global R Tc Pc Omega a b Kappa1 Kappa0 
R = 0.00008314;
Tc = 670.2; %in K
Pc = 55.02; %in Bar
Tboil = 434.7;
T = linspace(236.65,670,100);
P = linspace(0.0001,54.8987,100);
% T = 150;
% P = 40;
Tr = T/Tc;
Omega = 0.39983;
Kappa1 = -0.03471;
Kappa0 = 0.378893 + 1.4897153*Omega - 0.17131848*Omega^2 + 0.0196554*Omega^3;

for m = 1:length(T)
    for j = 1:length(P)
        Kappa = Kappa0 + Kappa1.*(1 + Tr.^0.5).*(0.7-Tr);
        DT1 = 1./T-1/Tboil;
        DT2 = 1/Tc-1/Tboil;
        DlogP = log(Pc);
        VaporPint = exp(DlogP.*DT1./DT2);
        b = 0.07780*R*Tc/Pc;
        AlphaSqrt = 1+Kappa.*(1-sqrt(Tr));
        Alpha = AlphaSqrt.^2;
        a = 0.45724*(R^2)*(Tc^2).*Alpha/Pc;
        Pv = (test5(T, VaporPint,a,b,R))  
    end
end
plot(T, VaporPint)

function of test5:

function Pv = test5(T,P,a,b,R)

criterion = 1;
K = 0;
while criterion >= 0.00001
    A = a.*P/((R.*T).^2);
    B = b*P/(R*T);
    Z = Zroots(A,B);
    Z(1) = max(Z);
    Z(2) = min(Z);
    fugacity = Solvefug(A,B,Z,P);
    criterion = abs(fugacity(2)/fugacity(1)-1);
    P = P*fugacity(2)/fugacity(1);
    K = K+1;
end
Pv = P;
fl = fugacity(2)
fv = fugacity(1)
end

My program works when I plug in one Temperature and one pressure, but when I try to plug in a range or T's and P's using linspace, I get an error. I have tried for hours to fix this issue but my efforts are sterile. I really need help. Any advice is helpful.

Thanks

[Walter Roberson]

Your test5 counts iterations in K but never uses K, so it does not seem useful to count them.

Answer 1

[Walter Roberson]

When you have

for m = 1 : 10
  a = m^2;    %for example
end

then in every iteration of the for loop, you overwrite the complete variable "a". After the loop "a" would have only the last value that was assigned to it, the same as if you had only done the loop once

for m = 10
  a = m.^2;
end

If you want to use the values after the loop, such as to plot them, then you need to store all of the values:

for m = 1 : 10
  a(m) = m^2;    %for example
end
plot(a)