function [ys,check,M_params] = model_steadystate(ys,exe)

global M_

NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end

check = 0;

GY = 0.18;
MY = 0.15*4;
BY = 0.10*4;
EMBI = 1.0^0.25;
TBY = eps;
OMG = 0.55;
OM = 0.2/0.55;

MZSS = 1.038^0.25;
MVSS = 1;

p = 1.0^0.25;
PSS = p;

ax=1;
an=1;
mun = 1;
mux = 1;
h = 0.2;

rpn = 1;
rpx = 1;
rpt = 1;
rpm = 1;
rpms = 1;

GYW = 0.20;
BYW = 0.45*4;
MYW = 0.56*4;
hs = 0.2;
OMW = 0.00001;


DZSS = MZSS*(MVSS^inv(1-THETA));
DZTSS = MZSS*(MVSS^(THETA/(1-THETA)));

YXY = OMG;
YNY = 1-YXY;

mv = MVSS;
mz = MZSS;
dz = DZSS;
dzt = DZTSS;
RR = (dzt^SIG)/(BETTA);

PXX = PSS;
PMM = PSS;

pn = PSS;
px = PSS;
pm = PSS;
pt = PSS;

r = p * RR;
RSS = r/EMBI;
rs = RSS;
BR = r;


q = inv(r);
qs = inv(rs);
PFSS = BETTA  / (qs * DZTSS^SIG);
pf = PFSS;
rf = r*pf/p;
epso = (rf/rs)^inv(K1);
EPSOSS = epso;
pxs = PFSS;
pms = PFSS;
PMSS = PSS;

mcw = GAMB/(GAMB-1);


pntil = ((1-ALFAN*pn^(ZETAN-1))/(1-ALFAN))^inv(1-ZETAN);
pxtil = ((1-ALFAN*px^(ZETAN-1))/(1-ALFAN))^inv(1-ZETAN);
pmtil = ((1-ALFAM*pm^(ZETAM-1))/(1-ALFAM))^inv(1-ZETAM);
pxstil = ((1-ALFAXP*pxs^(ZETAXP-1))/(1-ALFAXP))^inv(1-ZETAXP);
ptilf = ((1-ALFAW*pf^(ZETAW-1))/(1-ALFAW))^inv(1-ZETAW);

SN = ((1-ALFAN)*pntil^(-ZETAN))/(1-ALFAN*pn^ZETAN);
SX = ((1-ALFAN)*pxtil^(-ZETAN))/(1-ALFAN*px^ZETAN);
SM = ((1-ALFAM)*pmtil^(-ZETAM))/(1-ALFAM*pm^ZETAM);
SXP = ((1-ALFAXP)*pxstil^(-ZETAXP))/(1-ALFAXP*pxs^ZETAXP);    

sn = SN;
sx = SX;
sm = SM;
sxp = SXP;
sf = ((1-ALFAW)*ptilf^(-ZETAW))/(1-ALFAW*pf^ZETAW);

mcx = (pxtil*(1-ALFAN*q*dzt*(px^(KX-(1+ZETAN)/ZETAN))^(-ZETAN)))*(ZETAN-1)/(ZETAN*(1-ALFAN*q*dzt*(px^(KX-ZETAN/(ZETAN-1)))^(1-ZETAN)));
mcn = (pntil*(1-ALFAN*q*dzt*(pn^(KN-(1+ZETAN)/ZETAN))^(-ZETAN)))*(ZETAN-1)/(ZETAN*(1-ALFAN*q*dzt*(pn^(KN-ZETAN/(ZETAN-1)))^(1-ZETAN)));
mcs = (ptilf*(1-ALFAW*qs*dzt*(pf^(KW-(1+ZETAW)/ZETAW))^(-ZETAW)))*(ZETAW-1)/(ZETAW*(1-ALFAW*qs*dzt*(pf^(KW-ZETAW/(ZETAW-1)))^(1-ZETAW)));
rer = (ZETAM-1)*pmtil*(rpm*(1-ALFAM*q*dzt*(pm^(KM-(1+ZETAM)/ZETAM))^(-ZETAM)))/((1-ALFAM*q*dzt*(pm^(KM-ZETAM/(ZETAM-1)))^(1-ZETAM))*ZETAM);
RERSS = rer;
tot = (rpx*rpt)*(ZETAXP*(1-ALFAXP*q*dzt*(pxs^(KXP-ZETAXP/(ZETAXP-1)))^(1-ZETAXP)))/(rer*pxstil*(ZETAXP-1)*(1-ALFAXP*q*dzt*(pxs^(KXP-(1+ZETAXP)/ZETAXP))^(-ZETAXP)));
TOTSS = tot;

LY = MY + r*BY;
taud = GY - (r-1)*MY/r + (r-p*dzt)*LY/(r*p*dzt);
TAUD = taud;

rkn = inv(1-taud)*(inv(BETTA)*mv*dzt^SIG - 1 + DELTA);
rkx = inv(1-taud)*(inv(BETTA)*mv*dzt^SIG - 1 + DELTA);
khx = (rkx/(mcx*THETA))^(1/(THETA-1))*dz;
HXN = (mcx*YXY)/(mcn*YNY);
hn = h*inv(1+HXN);
hx = HXN*hn;

kbx = khx * hx;
kbn = (kbx*hn/hx) * (mcx/mcn)^(1/THETA);
w = mcx*(1-THETA)*((kbx/(hx*dz))^THETA);

ixd = (1-(1-DELTA)/dz)*kbx;
ivnd = (1-(1-DELTA)/dz)*kbn;
ive = ixd + ivnd;

T1 = (1-taud)*rkx;
T2 = T1*TK;

g = GY*(w*h + rkn*kbn/dz + rkx*kbx/dz);
GSS = g;
gt = OMG * g;
gn = (1-OMG) * g;
li = LY*(w*h + rkn*kbn/dz + rkx*kbx/dz);
LI = li;

c = (1-TBY)*(w*h + rkn*kbn/dz + rkx*kbx/dz)-ive-g;
cn = (1-OMG)*c;
ct = OMG*c;
cx = (1-OM)*ct;
cm = OM*ct;
ivn = (1-OMG)*ive;
it = OMG*ive;
ix = (1-OM)*it;
im = OM*it;
dmm = sm*(im + cm);
vx = 1;
vn = 1;

m = MY*(w*h + rkn*kbn/dz + rkx*kbx/dz);
VM = m/c;
IBY = inv(rer*(rf/dzt - pf))*TBY;
ibx = IBY*(w*h + rkn*kbn/dz + rkx*kbx/dz);
xs = (TBY * (w*h + rkn*kbn/dz + rkx*kbx/dz) + dmm);
dxp = sxp*xs;
IB = ibx;
cw = (xs*tot^ETW)/OMW;

vw = 1;
LYW = MYW + rs*BYW;
taus = (GYW - (rs-1)*MYW/rs + (rs-pf*dzt)*LYW/(rs*pf*dzt));
rkw = inv(1-taus)*(inv(BETTA)*mv*dzt^SIG - 1 + DELTA);
kw = (rkw/(mcs*THETA))^(1/(THETA-1))*hs*dz;
ws = mcs*(1-THETA)*((kw/(hs*dz))^THETA);
ivw = (1-(1-DELTA)/dz)*kw;
yw = (cw + cm + im + ivw)/(1 - GYW);
gw = GYW*yw;
GW = gw;
mf = MYW*yw;
MFSS = mf;
lis = LYW*yw;
bs = BYW*yw;
TAUS = taus;
VMW = mf/cw;
KKF = kw;
YHW = yw;
FCF = ((kw/dz)^THETA)*(hs^(1-THETA)) - sf*(yw); 
GAMW = (ws * (1 - taus) * (1-hs)) / ((cw-QSIW*cw/dzt) * (1 + VMW*(rs - 1)/rs));
lambdw = ((cw-QSIW*cw/dzt)^(-SIG) * (1 - hs)^(GAMW*(1 - SIG)))/(1 + VMW*(rs - 1)/rs); 


FCX = ((kbx/dz)^THETA)*(hx^(1-THETA))-sx*(cx+gt+ix+dxp);
FCN = ((kbn/dz)^THETA)*(hn^(1-THETA))-sn*(cn+gn+ivn);
FCM = dmm-sm*(cm+im);
FCXP = dxp - sxp*xs;

GAM = (w * (1 - taud) * (1-h)) / ((c-QSI*c/dzt) * (mcw*(1 + VM*(r - 1)/r)));
lambd = ((c-QSI*c/dzt)^(-SIG) * (1 - h)^(GAM*(1 - SIG)))/(1 + VM*(r - 1)/r); 

zt= (cx+gt+ix+dxp)*(ZETAN-1)*(pxtil^(-ZETAN))/((1-ALFAN*q*dzt*(px^(KX-ZETAN/(ZETAN-1)))^(1-ZETAN))*ZETAN);
zo= (cx+gt+ix+dxp)*mcx*(pxtil^(-1-ZETAN))/(1-ALFAN*dzt*q*(px^(KX-(1+ZETAN)/ZETAN))^(-ZETAN));
xt= (cn+gn+ivn)*(ZETAN-1)*(pntil^(-ZETAN))/((1-ALFAN*q*dzt*(pn^(KN-ZETAN/(ZETAN-1)))^(1-ZETAN))*ZETAN);
xo= (cn+gn+ivn)*mcn*(pntil^(-1-ZETAN))/(1-ALFAN*q*dzt*(pn^(KN-(1+ZETAN)/ZETAN))^(-ZETAN));
yt= (cm+im)*(ZETAM-1)*(pmtil^(-ZETAM))/((1-ALFAM*q*dzt*(pm^(KM-ZETAM/(ZETAM-1)))^(1-ZETAM))*ZETAM);
yo= (cm+im)*rer*(pmtil^(-1-ZETAM))/(1-ALFAM*q*dzt*(pm^(KM-(1+ZETAM)/ZETAM))^(-ZETAM));
uo= (pxstil^(-1-ZETAXP))*xs*(rpx*rpt)/(rer*tot*(1-ALFAXP*q*dzt*(pxs^(KXP-(1+ZETAXP)/ZETAXP))^(-ZETAXP)));
ut= (pxstil^(-ZETAXP))*xs*(ZETAXP-1)/(ZETAXP*(1-ALFAXP*q*dzt*(pxs^(KXP-ZETAXP/(ZETAXP-1)))^(1-ZETAXP)));
xwt = (cw+gw+cm+im+ivw)*(ZETAW-1)*(ptilf^(-ZETAW))/((1-ALFAW*qs*dzt*(pf^(KW-ZETAW/(ZETAW-1)))^(1-ZETAW))*ZETAW);
xwo = (cw+gw+cm+im+ivw)*mcs*(ptilf^(-1-ZETAW))/(1-ALFAW*dzt*qs*(pf^(KW-(1+ZETAW)/ZETAW))^(-ZETAW));


y = c + ive + g + rpx*rpt*xs - (rpm*rpt)*(dmm);
YH = y;
YH_CE = y;

Util = (((c - QSI*c/dzt)*((1-h)^GAM))^(1-SIG)-1)/(1-SIG);
Welf = Util/(1-BETTA); 
etam = ZETAM;
etax = ZETAN;
etan = ZETAN;
etaxp = ZETAXP;
etaw = ZETAW;


sx =log(sx);
c = log(c);
h = log(h);
w  = log(w);
p = log(p);
rpn = log(rpn);
rpx = log(rpx);
r = log(r);
mcn = log(mcn);
mcx = log(mcx);

xo = log(xo);
xt = log(xt);
q = log(q);
pxtil = log(pxtil);
ax = log(ax);
an = log(an);
cx = log(cx);
cn = log(cn);
px = log(px);
pn = log(pn);

hx = log(hx);
hn = log(hn);
zo = log(zo);
zt = log(zt);
pntil = log(pntil);
sn = log(sn);
lambd = log(lambd);
taud = log(taud);
g = log(g);
gt = log(gt);

gn = log(gn);
rkn = log(rkn);
rkx = log(rkx);
kbn = log(kbn);
kbx = log(kbx);
ix = log(ix);
ivn = log(ivn);
vx = log(vx);
vn = log(vn);
ive = log(ive);

ixd = log(ixd);
ivnd = log(ivnd);
mux = log(mux);
mun = log(mun);
m = log(m);
yw = log(yw);
pf = log(pf);
rf = log(rf);
epso = log(epso);
rs = log(rs);

rer = log(rer);
xs = log(xs);
ibx = log(ibx);
mv = log(mv);
mz = log(mz);
dz = log(dz);
dzt = log(dzt);
sm = log(sm);
im = log(im);
it = log(it);

cm = log(cm);
ct = log(ct);
pmtil = log(pmtil);
yo = log(yo);
yt = log(yt);
rpt = log(rpt);
rpm = log(rpm);
pm = log(pm);
pt = log(pt);
y = log(y);

pxs = log(pxs);
tot = log(tot);
pxstil = log(pxstil);
sxp = log(sxp);
dxp = log(dxp);
uo = log(uo);
ut = log(ut);
rpms = log(rpms);
pms = log(pms);
mcw = log(mcw);

mf = log(mf);
dmm = log(dmm);
etam = log(etam);
etax = log(etax);
etan = log(etan);
etaxp = log(etaxp);
etaw = log(etaw);
gw = log(gw);
lis = log(lis);
li = log(li);

taus = log(taus);
ws = log(ws);
hs = log(hs);
mcs = log(mcs);
sf = log(sf);
ptilf = log(ptilf);
xwt = log(xwt);
xwo = log(xwo);
cw = log(cw);
qs = log(qs);
bs = log(bs);
kw = log(kw);
ivw = log(ivw);
rkw = log(rkw);
vw = log(vw);

lambdw = log(lambdw);

NumberOfParameters = M_.param_nbr;                             % Number of parameters.
M_.params = zeros(NumberOfParameters,1);                       % Initialization of vector of parameters.
for i = 1:NumberOfParameters                                   % Loop...
  varname = deblank(M_.param_names(i,:));                      %    Get the name of parameter i.                     
  eval(['M_.params(' int2str(i) ') = ' varname ';']);          %    Get the steady state value of this variable.
end

NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for i = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
  eval(['ys(' int2str(i) ') = ' varname ';']);                %    Get the steady state value of this variable.
end