    
function F = findSS(X)
                
% Function to compute steady state for the model with unseparable utility.
% This is done numerically because of most of the variables are determined
% simul.
    
%X(1) = QHTOC;
%X(2) = QH1TOC1;
%X(3) = n;
%X(4) = n1;
%X(5) = CTOY;
%X(6) = C1TOY;
%X(7) = y;
%X(8) = q;
%X(9) = h;
%X(10) = h1;

paramfile_baby00

r = PIBAR / BETA;
rk = PIBAR / BETA - (1-DK);
xp = XP_SS;
xw = XW_SS;
xw1 = XW_SS;
lm = (1-BETA1/BETA) / (1-((BETA1*RHOB)/PIBAR));

KTOY = ALPHA / (xp*rk); 
B1TOQH1 = (M*(1 - RHOB)) /(1-(RHOB/PIBAR)); 
WNTOY = (1-ALPHA)*(1-SIGMA)/xp;
TTOY = (1/BETA-1)*BGTOY+GTOY;


F = [ -X(1)+((1/RHOD)*(1/((1-X(3))^((1-RHOD)*(1-VARTHETA)))*(JEI/(1-BETA))))^(-1/(RHOD*(1-VARTHETA)-1))*X(8)^(-RHOD*(1-VARTHETA)-1);
      -X(2)+((1/RHOD)*(1/((1-X(4))^((1-RHOD)*(1-VARTHETA)))*(JEI/(1-BETA1-lm*(1-RHOB)*M))))^(-1/(RHOD*(1-VARTHETA)-1))*X(8)^(-RHOD*(1-VARTHETA)-1);
      -X(3)+1/(((1+(xp*xw*((1-RHOD)/RHOD)*X(5))/((1-SIGMA)*(1-ALPHA)))));
      -X(4)+1/(((1+(xp*xw1*((1-RHOD)/RHOD)*X(6))/((1-SIGMA)*(1-ALPHA)))));
      -X(5)+(1-X(6)-DK*KTOY-GTOY);  
      -X(6)+((1-ALPHA)*SIGMA*(1/xp) - SIGMA*TTOY)/(1+(1/BETA-1)*B1TOQH1*X(2));
      -X(7)+KTOY^(ALPHA/(1-ALPHA))*X(3)^(1-SIGMA)*X(4)^(SIGMA);
      -X(8)+X(1)*X(5)*X(7)+X(2)*X(6)*X(7);
      -X(9)+(X(1)*X(5)*X(7))/X(8);
      -X(10)+(X(2)*X(6)*X(7))/X(8);
      ];
end

 