function [ys,check] = model1_steadystate(ys,exo)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%% Start SS equations here

r=1/beta_s;

muh=1-beta_h*r;
QHh_Ch=1/(1-beta_h-theta_h*muh);
WLh_Ch=1+(r-1)*theta_h*QHh_Ch;
psi_h=WLh_Ch;
lah=(WLh_Ch/psi_h)^(1/(1+etta));

mul=1-beta_i*r;
Kl_PlYl=alfa*(1-gama)*beta_i/(1-beta_i*(1-deltta)-theta_l*mul);
QHl_PlYl=alfa*gama*beta_i/(1-beta_i-theta_l*mul);

mue=1-beta_i*r;
Ke_PeYe=alfa*beta_i*(1+theta_e*mue)/(1-beta_i*(1-deltta));

F1 = @(x)[(lah+x(1))*x(2)-(1-alfa)/(psi_h*lah^etta);
          (1+psi_h*lah^etta/(x(3)*x(1)^etta))*x(2)-1+deltta*(chhi*Kl_PlYl+(1-chhi)*Ke_PeYe);
          x(1)-1];
x0=[1,0.190471566853713,0.485253226606126];
options = optimoptions('fsolve','Display','off','MaxFunctionEvaluations',10000,'MaxIterations',10000,'OptimalityTolerance',1e-25,'FunctionTolerance',1e-25,'StepTolerance',1e-25);
[ss1,fval1,exitf1]=fsolve(F1,x0,options);
if exitf1<1
    check=1;
end

las=ss1(1);
Ch_Y=ss1(2);
psi_s=ss1(3);

QHh_Y=QHh_Ch*Ch_Y;
QHl_Y=chhi*QHl_PlYl;
Hh_Hl=QHh_Y/QHl_Y;
hl=HH/(1+Hh_Hl);
hh=HH-hl;

Ll_Le=chhi/(1-chhi);
le=(las+lah)/(1+Ll_Le);
ll=las+lah-le;

XX=le^alfa*hl^(alfa*gama)*(Kl_PlYl*hl^(alfa*gama)*ll^(1-alfa))^(alfa*(1-gama)/(1-alfa*(1-gama)));
YY=ll^alfa*(Ke_PeYe*le^(1-alfa))^(alfa/(1-alfa));

F2 = @(y)[y(2)^(1/(1-alfa))/y(1)^(1/(1-alfa*(1-gama)))-XX/YY;
          chhi^chhi*(1-chhi)^(1-chhi)-y(1)^chhi*y(2)^(1-chhi)];
y0=[0.456298078191519,0.547887470819172];
[ss2,fval2,exitf2]=fsolve(F2,y0,options);
if exitf2<1
    check=1;
end

pl=ss2(1);
pe=ss2(2);

kl=(pl*Kl_PlYl*hl^(alfa*gama)*ll^(1-alfa))^(1/(1-alfa*(1-gama)));
ke=(pe*Ke_PeYe*le^(1-alfa))^(1/(1-alfa));

yl=(hl^gama*kl^(1-gama))^alfa*ll^(1-alfa);
ye=ke^alfa*le^(1-alfa);
y=yl^chhi*ye^(1-chhi);

ch=Ch_Y*y;
cs=y-ch-deltta*(ke+kl);

w=psi_h*lah^etta*ch;
q=ch/(hh*(1-beta_h-theta_h*muh));
pie=pe*ye-w*le;

bh=theta_h*q*hh;
bl=theta_l*(q*hl+kl);
be=theta_e*pie;
b=bh+bl+be;

cs=log(cs);
ch=log(ch); 
las=log(las);
lah=log(lah);
ll=log(ll);
le=log(le);
hh=log(hh);
hl=log(hl);    
bh=log(bh);
bl=log(bl);
be=log(be);
b=log(b);
yl=log(yl);
ye=log(ye);
y=log(y);
il=log(deltta*kl);
ie=log(deltta*ke);
kl=log(kl);
ke=log(ke);
pie=log(pie); 
w=log(w);
r=r;
q=log(q);
qkl=0;
qke=0;
pl=log(pl);
pe=log(pe);
mgus=-cs;
mguh=-ch;    
lams=0;
lamh=0;
muh=log(muh);
mul=log(mul);
mue=log(mue);
    
ea=0;
eh=0;
ei=0;


% End SS equations here

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end

end

