function [ys,params,check] = Model_steadystate(ys,exo,M_,options_)

check = 0;
options = optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

for ii = 1:M_.param_nbr
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
end

e0=4.02;
e=dsge_steadystate_helper(e0,rhoe,epsilone);
a0=1;
a=dsge_steadystate_helpera(a0,rhoa,epsilona);
x0=1;
x=dsge_steadystate_helperx(x0,epsilonx,rhox);
z0=3394.8;
z=dsge_steadystate_helperz(z0,rhoz,epsilonz);
v0=1;
v=dsge_steadystate_helperv(v0,epsilonv,rhov);
q=1/bet-1+delt;
lambd=(et+(1-alph)*((1+e*(r/(r-1))^(gamm-1))^(-1))*(thet/(thet-1)-delt*(alph/q))^(-1))/((1-alph)*z*((thet-1)/thet)^((1+alph)/(1-alph))*(alph/q)^(alph/(1-alph)));
c=((1+e*((r/(r-1))^(gamm-1))^(-1))*(1/lambd));
m=c*((r/(r-1))^(gamm))*c;
tau=dsge_steadystate_helpert(tau0,m,pi);
pi=tau;
r=pi/bet;
ksi=((thet-1)/thet)*lambd;
y=(1-delt*(alph/q)*((thet-1))/thet)^(-1);
k=(alph*q)*((thet-1)/thet)*y;
i=delt*k;
h=(1/z)*(y/k^alph)^(1/(1-alph));
w=(1-alph)*((thet-1)/thet)*(y/h);
wd=y-w*h-q*k;
n=y/h;

params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

for ii = 1:M_.orig_endo_nbr
  eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
end