% Find the steady state for chapter2_mod
function xx = chpater2_ss(xx)

global An phi0 sigma r am; 
global alphan alpha nu phi Ome epsilon etay etae eta lambda rho;
global pmk_ss pmc_ss px_ss pn_ss;
global lh;


px = px_ss;
pmk = pmk_ss;
pmc = pmc_ss;
pn = pn_ss;

  
zave = (px*etay+pn*eta)*(phi/(epsilon-phi));
xi = etae*(r+lambda)*pn/zave;
varphistar = phi0*(xi^(-1/epsilon));
varphiave = varphistar*((epsilon/(epsilon-phi))^(1/phi));
wx = (px*varphistar*((Ome/(px*etay+pn*eta))^(1/phi))/(pmk^nu))^(1/alpha);


%Based on wx = wn
delta = ((((px/(pmk^nu))*phi0*(Ome^(1/phi))* ...
            ((pn*etae*(r+lambda))^(-1/epsilon))* ...
            ((phi/(epsilon-phi))^(1/epsilon))* ...  %%% the part that missing
            (((px*etay)+(pn*eta))^((phi-epsilon)/(phi*epsilon)))/ ...
            ((pn*alphan*(An^(1/alphan)))^(alpha)))^(alphan/(alpha*(alphan-1))))*(1-alphan)) - r ;


rn = pn*(r+delta);


% cn = yn - in - etae*he - xi*eta*h
% cn = cnln*ln - cn*h from market clear condition of non-tradable goods
cnln = (An*(((1-alphan)*wx/(alphan*rn))^(1-alphan))) ...
           -(delta*(1-alphan)*wx/(alphan*rn));
cnh  = (etae*lambda) + (xi*eta);

% cm = (px*yx - pmk*kx)/pmc  from trade balance = 0
% thus cm = (wx*lx + xi*zave*h + xi*eta*h)/pmc
% since xi*zave*h = px*yx - pmk*kx - wx*lx - xi*eta*h
% cm = cmh*h
cmh = ((wx*(((px*varphiave)^phi)*((alpha/wx)^((1-nu)*phi))*...
          ((nu/pmk)^(nu*phi))))+zave+eta)*xi/pmc;       %%% /pmc

% ln = lntoh * h 
lntoh = (cnh+(((1-am)/am)*((pn/pmc)^(-rho))*cmh))/cnln;

%lx = lxtoh *h
lxtoh = xi*(((px*varphiave)^phi)*...
             ((alpha/wx)^((1-nu)*phi))*((nu/pmk)^(nu*phi)));

% l = ln + lx = ltoh * h
ltoh = lntoh+lxtoh;

% share of labor in two sectors
shln = lntoh/ltoh;
shlx = lxtoh/ltoh;

l = lh;
h = l/ltoh;

ln = lntoh*h;
kn = ln*(1-alphan)*wx/(alphan*rn);
yn = An*(ln^alphan)*(kn^(1-alphan));
in = delta*kn;

lx = lxtoh*h;
kx = lx*wx*nu/(alpha*pmk);
yx = ((wx*lx) + (pmk*kx) + (xi*zave*h) + (xi*pn*eta*h) + (px*etay*xi*h))/px;
he = lambda*h;

cn = yn-in-(etae*he)-(xi*eta*h);
cm = (((px*yx)-(px*etay*xi*h))-(pmk*kx))/pmc;

c = ((((1-am)^(1/rho))*(cn^((rho-1)/rho)))+...
      ((am^(1/rho))*(cm^((rho-1)/rho))))^(rho/(rho-1));
p = (((1-am)*(pn^(1-rho)))+(am*(pmc^(1-rho))))^(1/(1-rho));
%% check for the utility  
%u=((((c^kappa)*((1-l)^(1-kappa)))^(1-sigma))-1)/(1-sigma)

kappa = 1/((wx*(1-l)/(p*c))+1);   %must smaller than 1

varsigma = r/((c^kappa)*((1-l)^(1-kappa)));



xx(1)=delta;
xx(2)=kappa;
xx(3)=varsigma;
xx(4)=shln;
xx(5)=varphiave;


gdp = (p*c)+(pn*(in+(etae*he)));
k = kn + h;