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 trying as follows. My concern is from "for loop o=1:N". Please i could get some suggestion to write it properly. Mathematically, we got to solve it like numerical methods but in writing code I am not able to.

%%To calculate system pressure and vapor phase composition
  %y=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));  
    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