How to solve non-linear equations with several unknowns and several equations?

I have a problem with xz. After I run this code I get the same value of "xz" for 500 times.
"xz" should not get the same value. It should vary with P, T_sat and q. 
I think the problem is "eps" in the second round cannot get beyond "eps>0.0001" then the rest of "xz" has the same values. 

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));

     for i=2:501   
         eps=1;
         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

[John D'Errico]

It is a really bad idea to name a variable eps. This is a useful function in MATLAB.

For this code though, only you know what it is doing. Learn to use the debugger.

[Krittiya Tagong]

To be honest, I have tried by myself but I am struggling to solve it. If you think "eps" is a bad idea, would you mind suggest me which one I should use with this problem?

Thank you