Iteration using for loop and if statement

I am writing a code where I need to calculate P( as P_Initial) with all the variables calculated in other functions, then using this P_initial calculte PhiN, assign PhiN to Phi and once again calculate P(as P_total) with new value of Phi. Now if condition needs to be applied if true then print output ,if false recalculte PhiN using recent P_total and then check for condition again. This is what I need to do but while writing it in loop and to assign new value iteratively I am getting some issues. I am tring as follows. My concern is from "for loop o=1:N". Please i could get some suggestion to write it properly

%%To calculate system pressure and vapor phase composition
%yi=vapor composition
%P=Pressure(From Antoine's Equation)
%e=tolerance

e=0.4
[f_sat]=FugacityS(T,P,Tc,Pc)
[Gamma_t]=Gamma_sum(taua,x1,x2,x3,rho,zc,za,xc,xa,xm,T)
[Phi]=Redlich_kwong(T,Tc,Pc)
N=length(T);
M=length(Pc);
xm=[0.5096 0.5092 0.5087
    0.0963 0.0964 0.0965 
    0.3941 0.3944 0.3948];
 a0=zeros(1,M);aHI=zeros(1,N);aH2O=zeros(1,N);aI2=zeros(1,N);a=zeros(N,1);
for i=1:M
  a0(i)=(0.42748*(((R^2)*(Tc(i)^2.5))/Pc(i)));
end
for j=1:N
  for i=1
    aHI(j)=(((a0(i)*a0(i+1))^0.5)*xm(i,j)*xm(i+1,j))+(((a0(i)*a0(i+2))^0.5)*xm(i,j)*xm(i+2,j));  %applying mixing rules
    aH2O(j)=(((a0(i+1)*a0(i))^0.5)*xm(i+1,j)*xm(i,j))+(((a0(i+1)*a0(i+2))^0.5)*xm(i+1,j)*xm(i+2,j));
    aI2(j)=(((a0(i+2)*a0(i))^0.5)*xm(i+2,j)*xm(i,j))+(((a0(i+2)*a0(i+1))^0.5)*xm(i+2,j)*xm(i+1,j));
  end
end
format long g
for j=1:N
     a(j)=(aHI(j)+aH2O(j)+aI2(j));
     a(j)=(a(j))';
end
b0=zeros(N,M);b=zeros(1,N);A0=zeros(N,M);A=zeros(1,N);B0=zeros(N,M);
B=zeros(1,N); 
for j=1:N
  for i=1:M
    b0(i,j)=(0.08664*((R*Tc(i))/Pc(i)))*(xm(i,j));
    A0(i,j)=(((0.4278*Tc(i)^2.5)/(Pc(i)*T(j)^2.5))^(1/2))*(xm(i,j));
    B0(i,j)=((0.08664*Tc(i))/(Pc(i)*T(j))*(xm(i,j)));
  end
  b(j)=sum(b0(:,j));
  A(j)=sum(A0(:,j));
  B(j)=sum(B0(:,j));
end
Asq=((A).^2);
P_calI=zeros(N,M); y=zeros(N,M); P_Initial=zeros(N,1);y_total=zeros(N,1);
for o=1:N
  for p=1:M    
    P_calI(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/Phi(p,o)
  end
  P_Initial(p,o)=sum(P_calI(:,o)) 
  Z=zeros(1,M);h=zeros(1,M);PhiN=zeros(1,M);   
  h(p,o)=b(o)./Vm(p,o);
  Z(p,o)=(1./1-h(p,o))-((Asq(o)/B(o)).*(h(p,o)./(1+h(p,o))));
  PhiN(p,o)=exp(Z(p,o)-1-log(Z(p,o)-(B(o).*P_cal(p,o)))-(Asq(o)/B(o)).*(log(1+((B(o).*P_calI(p,o))./Z(p,o)))))
  Phi=PhiN; P_total=zeros(N,M);P_cal=zeros(N,M)
  P_cal(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/Phi(p,o)
  P_total(p,o)=sum(P_cal(:,o)) 
  if ((ptotal-p_Initial/p_Initial)<e)
    y(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/(Phi(p,o)*P_Initial(p,o))
    disp(p_total)
    disp(y)
  else
    Z=zeros(1,M);h=zeros(1,M);PhiN=zeros(1,M);   
    h(p,o)=b(o)./Vm(p,o);
    Z(p,o)=(1./1-h(p,o))-((Asq(o)/B(o)).*(h(p,o)./(1+h(p,o))));
    PhiN(p,o)=exp(Z(p,o)-1-log(Z(p,o)-(B(o).*P_cal(p,o)))-(Asq(o)/B(o)).*(log(1+((B(o).*P_cal(p,o))./Z(p,o)))))
    Phi=PhiN; P_total1=zeros(N,M)
    P_cal(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/Phi(p,o)
    P_total1(p,o)=sum(P_cal(:,o)) 
   end
 end

[Bob Thompson]

"but while writing it in a loop and to assign new value iteratively I am getting some issues.'

What issues are you experiencing? We cannot run your code directly as posted because you have several undefined variables and functions, but we might be able to help you find the solution if you give us some more detail about what problems you're seeing.

[GAGANDEEP KAUR]

I understand your concern and I was also considering it quite odd but all other variables are calculated from 6-7 separate function files and needs to be called here which will be quite complicated to share here.

Still if any general commands can be suggested for writing, as we generally perform these steps in numerical methods,like to get initial guess( here P_initial) then from that PhiN and for second value of P(P_totoal here) recently calculated Phi(PhiN) then check error if lesser than accepted end by displying calculated variables if not calculate Phi again using recent P. Mthematically I have idea but in writing I am not able .

[GAGANDEEP KAUR]

Here, I am sharing quite edited version of aforementioned code with values of variables needed. Actually I am getting problem in re-entering the for loop need to enter p loop again where ym is to calculted again with new parameters calculated until condition is not satisfied.After the condition get satisfied, then it should enter O loop(of ym) and same calculations to be done.

%%To calculate system pressure and vapor phase composition

e=0.3

f_sat= [ 32.2587857596876 34.0501961493468 35.8746854026426

2.00528670528962 2.33558884889119 2.70783015514365

0.157445079883318 0.187083254240374 0.221158510598425];

Gamma_t=[3.19702443124609 2.54070700562388 2.08314683871606;

3.31850686110491 2.66546741032482 2.20838701999553;

3.44415347633114 2.79709939502664 2.34012176373733];

Vm= [ 68.0045214121783 70.0736822139801 72.5465660407398

18.1158800129555 18.2324752448415 18.3515283669796

51.5543766322469 51.7464896478157 51.9409677300663];

Pc= [831,220.58,116.54];

Tc=[394.15 399.15 404.15];

xm=[0.09630,0.0964,0.0965;

0.5096,0.5092,0.5087;

0.3941,0.3944,0.3948];

N=length(xm);

M=length(Pc);

Phi(1:3,1:3)=1;

P_calI=zeros(N,M);ym=zeros(N,M);P_0=zeros(N,1);

for o=1:N

for p=1:M

P_calI(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/Phi(p,o)

end

P_0(o)=sum(P_calI(:,o))

end

for o=1:N

for p=1:M

ym(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/(Phi(p,o)*P_0(o))

end

a0=zeros(1,M);aHI=zeros(1,N);aH2O=zeros(1,N);aI2=zeros(1,N);a=zeros(N,1);

R=0.08314;

for i=1:M

a0(i)=(0.42748*(((R^2)*(Tc(i)^2.5))/Pc(i)));

end

for i=1

aHI(o)=(((a0(i)*a0(i+1))^0.5)*ym(i,o)*ym(i+1,o))+(((a0(i)*a0(i+2))^0.5)*ym(i,o)*ym(i+2,o)); %applying mixing rules

aH2O(o)=(((a0(i+1)*a0(i))^0.5)*ym(i+1,o)*ym(i,o))+(((a0(i+1)*a0(i+2))^0.5)*ym(i+1,o)*ym(i+2,o));

aI2(o)=(((a0(i+2)*a0(i))^0.5)*ym(i+2,o)*ym(i,o))+(((a0(i+2)*a0(i+1))^0.5)*ym(i+2,o)*ym(i+1,o));

end

format long g

a(o)=(aHI(o)+aH2O(o)+aI2(o));

a(j)=(a(j))';

b0=zeros(N,M);b=zeros(1,N);A0=zeros(N,M);A=zeros(1,N);B0=zeros(N,M);

B=zeros(1,N);

for i=1:M

b0(i,o)=(0.08664*((R*Tc(i))/Pc(i)))*(ym(i,o));

A0(i,o)=(((0.4278*Tc(i)^2.5)/(Pc(i)*T(o)^2.5))^(1/2))*(ym(i,o));

B0(i,o)=((0.08664*Tc(i))/(Pc(i)*T(o))*(ym(i,o)));

end

b(o)=sum(b0(:,o));

A(o)=sum(A0(:,o));

B(o)=sum(B0(:,o));

Asq=((A).^2);

Z=zeros(N,M);h=zeros(N,M);Phi=zeros(N,M); P_1=zeros(N,1);P_cal=zeros(N,M)

for p=1:M

h(p,o)=b(o)./Vm(p,o)

Z(p,o)=(1./1-h(p,o))-((Asq(o)/B(o)).*(h(p,o)./(1+h(p,o))))

Phi(p,o)=exp(Z(p,o)-1-log(Z(p,o)-(B(o).*P_calI(p,o)))-(Asq(o)/B(o)).*(log(1+((B(o).*P_calI(p,o))./Z(p,o)))))

P_cal(p,o)= (xm(p,o)*Gamma_t(p,o)*f_sat(p,o))/Phi(p,o)

end

P_1(o)=sum(P_cal(:,o))

if ((P_1(o)-P_0(o)/P_0(o))<e)

break

end

P_0=P_1

P_calI=P_cal

end