How can I make a loop with multivariable within a while-loop

I tried to get xz which has 500 values. xz is related to q. q is related to Tsat. Tsat is related to P. and function of P is related to xz again. This loop works well with the first value, but there is a problem with the second value which could not achieve "eps>0.0001". I have no idea what I should do with this code. One more question, why the first value is zero ?

T0=-90;
%Tz=1*ones(500,1);
A0 = 5*10^-6;
z= -500:0;
gamma_0 = 0.4 ;
D =1200;
q0= 150*10^-3;
a= 0.001;
m_fluid = 150;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
m_steam = 30; 
d = 0.4318;
rho_l= 1000;
rho_g= 0.804;
G= (4*(m_fluid+m_steam))/(pi()*d^2);
u = (G/rho_l);
g=9.81;
alpha0=0.1;
xz(1) = rho_g*alpha0/((rho_l*(1-alpha0))+(rho_g*alpha0));
eps=1;

   for i=2:501   

       Tz(i) = -(1/a)*((1+a*T0)*exp(((a*A0*D.^2/gamma_0)*(1-exp(-z(i)/D))+((a*q0*z(i))/(gamma_0+gamma_0*a*T0)))-((a*A0*D*z(i))/gamma_0))-1);
      sigma(i)=((-7.5*(10^(-12)))*Tz(i).^4)+(5.8*(10^(-9))*(Tz(i).^3))-1.8*(10^(-6))*(Tz(i).^2)-(1.9*(10^(-5))*Tz(i))+0.07;
       h_lg(i) = ((-5.42*(10^(-5)))*Tz(i).^4)+((0.013)*(Tz(i).^3))-4.5*(Tz(i).^2)-(1.9*(10^(3))*Tz(i))+2.5*(10^(6));
       Cp(i)= (0.0001*Tz(i).^3)-(0.027*Tz(i).^2)+ 3.9*Tz(i)+(4*10.^3);
       k(i)=((-1.3853*(10^(-10)))*Tz(i).^4)+(8.7543*(10^(-8)))*(Tz(i).^3)-(2.43*(10^(-5))*(Tz(i).^2)-(0.0030633)*Tz(i))+0.54455;
       Tz_k(i)=Tz(i)+273.15;
      mu_l(i)= 2.414*10.^-5*10.^(247.8/(Tz_k(i)-140));
      mu_g(i) = (18.27*10.^-6)*((291.15+120)/(Tz_k(i)+120))*((Tz_k(i)/291.15).^(3/2));
      Pr(i)=  ((Cp(i)*mu_l(i))/k(i));

       Tz_new=Tz(i);
       sigma_new=sigma(i);
       h_lg_new=h_lg(i);
      Pr_new=Pr(i);
      k_new=k(i);
      mu_l_new=mu_l(i);
      mu_g_new=mu_g(i);
      xz_new=xz(i-1);

      q_new = 150*10^-3;%Initial condition
      h_new = 1*10^5; %initial condition    it's here$ the inital condition is wrong

   while eps>0.0001

             xz_ref= xz_new;
             sigma_ref=sigma_new;
             h_lg_ref= h_lg_new;
             Tz_ref=Tz_new;
             Pr_ref=Pr_new;
             mu_l_ref=mu_l_new;
             mu_g_ref=mu_g_new;
             k_ref=k_new;
             q_ref=q_new;
             h_ref=h_new;% it copy here

          %%Pressure%%
          %%Firction Pressure%%
          rho_x = ((xz_ref/rho_g)+((1-xz_ref)/rho_l)).^-1;
          WE = ((rho_l^2)*(u^2)*d)/(g*sigma_ref*rho_x);
          Fr=((rho_l^2)*(u^2))/(g*(rho_x.^2));
          H=((rho_l/rho_g).^0.91)*((mu_g_ref/mu_l_ref).^0.19)*((1-(rho_g/rho_l)).^0.7);
          F=(xz_ref.^0.78)*((1-(xz_ref.^2)).^0.24);
          Re_l =((G*d)/mu_l_ref);
          Re_g =((G*d)/mu_g_ref);
          fl =0.079/(Re_l.^0.25);
          fg = 0.079/(Re_g.^0.25);
          E=(1-(xz_ref.^2))+((xz_ref.^2)*((rho_l*fg)/(rho_g*fl)));  
          frict2 = E+((3.24*F*H)/((Fr.^0.045)*(WE.^0.035)));
          PL = 4*fl*(-z(i)/d)*(G.^2)*(1/(2*rho_l));
          P_frict = PL*frict2;

          alpha = (1+0.28*((1+xz_ref)/xz_ref)^0.64*(rho_g/rho_l)*(mu_l_ref/mu_g_ref).^0.07)^(-1);
          rho_h = (rho_l*(1-alpha))+(rho_g*alpha);
          P_static = rho_h*g;

          P_mom = (G^2)*((((1-xz_ref).^2)/(rho_l*(1-alpha)))+ ((xz_ref).^2)/(rho_g*alpha))-(G^2)*((((1-xz(i-1)).^2)/(rho_l*(1-alpha)))+ ((xz(i-1).^2)/(rho_g*alpha)));

          P= P_frict+P_static-P_mom;
          Tsat = -1.8*10.^(-24)*P.^4+2*10.^(-17)*P.^3-8.1*10.^(-11)*P.^2+0.00016*P+88;
          h_l = 0.023*(k_ref/d)*((G*(1-xz_ref)*d)/(mu_l_ref).^(0.8))*(Pr_ref.^(1/3));  
          q_new = h_ref*(Tz_ref-Tsat);
          h_new = h_l*((1+(3000*((q_ref/G*h_lg_ref).^0.86)+(((xz_ref/(1-xz_ref)).^(3/4))*((rho_l/rho_g).^0.41)))));      
          xz_new= xz_ref+((4*q_ref)/(d*G*h_lg_ref))*(-z(i-1)+z(i));
          eps = abs(xz_ref - xz_new);

        end
        xz(i) = xz_new;
        Tz(i) = Tz_new;
        sigma(i)=sigma_new;
        %Tz_k=Tz_k_new;
        mu_l(i)= mu_l_new;
        mu_g(i) = mu_g_new;
        h_lg(i)= h_lg_new;
        Pr(i)=Pr_new;
        q=q_new;
        h=h_new;

   end

Thank you

[Sanne]

just a thought spin: If xz is related to q, then why do you keep iterating over q0, a constant? I guess this constant value will always give you constant values for xz over i. The only thing that is changing per i is that the new value for xz is the old value for xz added to a function dependent on q (constant) and z(-500:0). Was that the intention?

[Krittiya Tagong]

q0 and q are different. q0 uses in Tz only and it's a constant. But for "q" is used for finding xz.

[Sanne]

Hmm, firstly when I run your code, it works for all i's. However, I get all the same values for xz (all close to zero and smaller than 0.001), which makes sense:

for i=...

xz_new=xz(i-1);

   while eps>0.0001
             xz_ref= xz_new;
             xz_new= xz_ref+((4*q_ref)/(d*G*h_lg_ref))*(-z(i-1)+z(i));
             eps = abs(xz_ref - xz_new);
   end

   xz(i) = xz_new;
end

I dont know what exactly youre trying to do, but maybe this is an ordinary differential equation problem?

As a side note, I would recommend preallocating your xz(i)'s, Tz(i)'s etc.

Best, Sanne

[Krittiya Tagong]

Yes, I got the same as you that "xz" has the same value for 500 times. This one is my problem too. In fact, this one should get 500 values of xz difference 500 times. Tz(i) works well which has 500 values and relates to z(-500,0).

For xz should iterate follow P, Tsat, and q. All the values are related together.