Seek Matlab Help in this problem

1. Feb 11, 2012

qazinasir

% Calculate the vapor Pressure of n-heptane using SRK-EOS
% P = RT / V - b - alpha*a / V(V-b)
R = 8.134e-3; % MPa m3 / mole.K
n = 0;
n = n+1;
for i=1:5
Pn = 1.6714; %(Assume Pressure)> this needs updating
Tr = 0.9;
Tc = 540.13; % K
Pc = 2.74; % Mpa;
T = Tr*Tc;
w = 0.350;
m = 0.480+1.574*w -0.176*w^2;
al=(1+m*(1-Tr^0.5))^2;
a = 0.42748*(((R*Tc)^2) /Pc);
b = 0.08664*(R*Tc / Pc);
% Z3 - Z2 + (A-B-B^2)Z-AB=0
A = (a*Pn*al)/(R*T)^2;
B = (b*Pn)/(R*T);
Zv=max(roots([1 -1 A-B-B^2 -A*B]));
Zl=min(roots([1 -1 A-B-B^2 -A*B]));
disp(Zv)
disp(Zl)
phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv )));
phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl )));
Fv=Pn*phiv;
Fl=Pn*phil;
disp(phiv)
disp(phil);
disp(Fv);
disp(Fl);
Error=Fl-Fv;
if (Error < 10^-4);
Pnew=Pn;
else
P(n+1) =Pn * (phil/phiv);
Pn = P(n+1);
end
end
Basically in this code I am trying to assume a pressure 'Pn' to calculate the results which satisfy this relation (Error<10^-4) if this satisfy's then our assume pressure is right else if not then we use the relation P(n+1) = Pn*(Phil/Phiv) to caluclate another pressure and update the above pressure (previously assume) with this one. I have error with this and need advice please help

2. Feb 11, 2012

Pythagorean

it's a very small chance people will read that much code and do that much thinking for you. Try to narrow down your specific theoretical or technical problem. Try debugging your code a little bit to chase the source down; you might even solve the problem if you spend some time trying to understand it before you ask.

3. Feb 11, 2012

qazinasir

Basically in this code I am trying to assume a pressure 'Pn' to calculate the results which satisfy this relation (Error<10^-4) if this satisfy's then our assume pressure is right else if not then we use the relation P(n+1) = Pn*(Phil/Phiv) to caluclate another pressure and update the above pressure (previously assume) with this one. but this cycle of Assuming and error checking is perform in a loop (if i am not wrong ) . How to update the previously Assume pressure with new one. thats the main question