function F = steadystate(x)
global beta beta2 gamma delta delta2 theta k nu fe fe2 fv fh fh2 tau tau2 zmin zmin2 ZD2 L L2 t1_bar t2_bar; 
global alpha eta eta2 RHO11 RHO12 RHO21 RHO22 SIGMA1 SIGMA2 SIGMA12;
global NVN NHN NHN2 ZVC_bar ZHC_bar TOL_bar CC2Q_bar Q_bar ddw_bar dvw_bar;
global dhw_bar ZHC2_bar RH_bar RH2_bar NN2_bar;
global RHS4 RHS5 ZD_bar ZV_bar ZH_bar ZH2_bar N_bar NV_bar NH_bar ND2_bar NH2_bar;

%t1_bar=1;
%t2_bar=1;
   
 F=[(1-beta*(1-delta))*fe/(beta*(1-delta))-(1-x(1)^(-k))*x(6)-(x(1)^(-k))*x(7)-(x(2)^(-k))*x(8);
    
    x(6)-(k/(k-theta+1))*fh*x(4)*(tau2^(theta-1))*(t2_bar^(theta))*(x(5)^(1-theta))*((x(1)/x(2))^(theta-1))*((x(1)^(k-theta+1))-1)/(x(1)^k-1);
   
    x(7)-(k/(k-theta+1))*fh*(x(3)^(1-theta))*x(4)*((x(1)/x(2))^(theta-1))*((tau2/tau)^(theta-1))*((t2_bar/t1_bar)^(theta))*(x(5)^(1-theta))+fv*x(3);  
   
    x(8)-((theta-1)/(k-theta+1))*fh;
       
    x(7)-((x(1)^k-1)/(x(1)^(k-theta+1)-1))*x(6)-((theta-1)/(k-theta+1))*fv*x(3); 
             
    (1-beta2*(1-delta2))*fe2/(beta2*(1-delta2))-(k/(k-theta+1))*fh2*(tau*x(5))^(theta-1)*(t1_bar^theta)*(x(9)/zmin2)^(1-theta)*(x(4)^(-1))-(x(9)^(-k))*((theta-1)/(k-theta+1))*fh2;
       
    x(4)-((x(2)*tau/(x(9)*tau2))^(theta-1))*((t1_bar/t2_bar)^theta)*(x(5)^(theta-1))*(x(3)^theta)*fh2/fh;
              
    (x(1)^(-k))*(((theta-1)*k*fh2/(k-theta+1))*(x(1)/x(9))^(theta-1)+fv)+(k*theta/(k-theta+1))*(x(9)^(-k))*(x(12)^(-1))*fh2-(k*theta/(k-theta+1))*(x(2)^(-k))*fh*(x(3)^(-1));    
        
    x(5)^(1-theta)-(x(3)^(1-theta)+((tau2*t2_bar)^(1-theta))*(x(2)^(theta-1-k))*x(12))/((1-x(1)^(-k))*(x(1)^(theta-1))*((x(1)^(k-theta+1)-1)/(x(1)^k-1))*x(12)+(x(1)^(theta-1-k))*x(12)*(tau*t1_bar*x(3))^(1-theta)+(x(9)^(theta-1-k))*(tau*t1_bar*x(3))^(1-theta));
        
    ((theta*k*fh/(k-theta+1))*t2_bar*(1-(t2_bar-1)*NH_bar*(x(10)^(1-theta))*x(5))*(x(10)^(theta-1)/x(3)))-1-((1-beta2)/(beta2*(1-delta2)))*fe2*(x(10)^(theta-1))/x(14);
    
    x(14)-(1/x(2))^(theta-1)*(tau2*t2_bar/x(3))^(theta-1)-((1/x(2))^k)*x(12);
    
    (theta*k*fh2/(k-theta+1))*t1_bar*(1-(t1_bar-1)*NH2_bar*(x(11)^(1-theta))-(t1_bar-1)*NV_bar*(x(11)*x(9)/x(1))^(1-theta))*(x(11)^(theta-1))*x(3)-1-((1-beta)/(beta*(1-delta)))*fe*(x(11)^(theta-1))/x(13);
    
    x(13)-(1-(1/x(1))^k)*(x(1)/x(9))^(theta-1)*((x(1)^(k-theta+1)-1)/(x(1)^k-1))*(tau*t1_bar*x(3))^(theta-1)-((1/x(1))^k)*(x(1)/x(9))^(theta-1)-(1/x(12))*((1/x(9))^k);
    
    (x(12)*(x(10)/x(11))^(theta-1))-x(14)/x(13)];

ZV_bar = x(1);
ZH_bar = x(2);
TOL_bar = x(3);
CC2Q_bar = x(4);
Q_bar = x(5);
ddw_bar = x(6);
dvw_bar = x(7);
dhw_bar = x(8);
ZH2_bar = x(9);
RH_bar = x(10);
RH2_bar = x(11);
NN2_bar = x(12);
RHS4 = x(13);
RHS5 = x(14);
