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

paramfile

nparams = size(M_.param_names,1);
for icount = 1:nparams
    eval(['params(icount,1) = ',M_.param_names{icount,:},';'])
end

check=0;
%% THIS BLOCK IS MODEL SPECIFIC.
%=========
%preambles
%=========
Xi1=xip1;Xi2=xip2;Tau1=taup1;Tau2=taup2;L1=lp1;L2=lp2;Z1=zp1;Z2=zp2;
Gamb1 = Tau1 + (1-btpar*(1-Xi1))*(1-Tau1)*alpn + (1-Tau1)*alpl;
Gamb2 = Tau2 + (1-btpar*(1-Xi2))*(1-Tau2)*alpn + (1-Tau2)*alpl;
Gamd1 = thetpt*Tau1 - thetpl*(1-btpar)*(1-Xi1)*(1-Tau1)*alpn;
Gamd2 = thetpt*Tau2 - thetpl*(1-btpar)*(1-Xi2)*(1-Tau2)*alpn;

options=optimset('Maxfuneval',10000,'Display','on');

%=========
%step1
%=========
Rb = 1/betah;
Rk = Rb - 1 + delk;

%=========
%step2
%=========
% mg1py1 mg2py1 mg3v1 mg4py1 bOpy1 kgOpy1 qhlOpy1
% xx(1)  xx(2)  xx(3) xx(4)  xx(5) xx(6)  xx(7)
ff=@(xx)([
    xx(5)*(Rb-1) + delg*xx(6) - Gamb1;
    Gamd1 + thetpl*xx(7) - xx(5);
    pip*(1-betag*omg2)/(1-omg2)/delg/xx(6) - xx(1) + xx(2);
    (1-betag*omg1)/(1-omg1) + Gamb1*xx(1) + Gamd1*xx(4) - (1-alpk-alpn)/alpg*(1-betag*(1-delg))/betag*xx(2)*xx(6) + (1-alpk-alpn)/alpg*varphig*xx(3);
    (1-betag*Rb)*xx(1) - xx(4);
    (1-betag*omg1)/(1-omg1) + Gamb1*xx(1) + Gamd1*xx(4) - (1-alpk-alpn)/alpn*btpar*xx(3) - (1-alpk-alpn)/alpn*thetpl*xx(7)*xx(4);
    (1-btpar)*xx(3) - thetpl*xx(7)*xx(4)
    ]);
xx=fsolve(ff,[1,1,1,1,1,1,1]);
mg1py1=xx(1); mg2py1=xx(2); mg3v1=xx(3); mg4py1=xx(4); bOpy1=xx(5); kgOpy1=xx(6); qhlOpy1=xx(7);
lfOlmlf1 = alpl/(1-alpk-alpn)/thetpl/qhlOpy1/mg4py1*((1-betag*omg1)/(1-omg1) + Gamb1*mg1py1 + Gamd1*mg4py1);
Lf1 = L1*lfOlmlf1/(1+lfOlmlf1);
qhOpy1 = qhlOpy1/(L1-Lf1);
% mg1py2 mg2py2 mg3v2 mg4py2 bOpy2 kgOpy2 qhlOpy2
% yy(1)  yy(2)  yy(3) yy(4)  yy(5) yy(6)  yy(7)
gg=@(yy)([
    yy(5)*(Rb-1) + delg*yy(6) - Gamb2;
    Gamd2 + thetpl*yy(7) - yy(5);
    pip*(1-betag*omg2)/(1-omg2)/delg/yy(6) - yy(1) + yy(2);
    (1-betag*omg1)/(1-omg1) + Gamb2*yy(1) + Gamd2*yy(4) - (1-alpk-alpn)/alpg*(1-betag*(1-delg))/betag*yy(2)*yy(6) + (1-alpk-alpn)/alpg*varphig*yy(3);
    (1-betag*Rb)*yy(1) - yy(4);
    (1-betag*omg1)/(1-omg1) + Gamb2*yy(1) + Gamd2*yy(4) - (1-alpk-alpn)/alpn*btpar*yy(3) - (1-alpk-alpn)/alpn*thetpl*yy(7)*yy(4);
    (1-btpar)*yy(3) - thetpl*yy(7)*yy(4)
    ]);
yy=fsolve(gg,[1,1,1,1,1,1,1]);
mg1py2=yy(1); mg2py2=yy(2); mg3v2=yy(3); mg4py2=yy(4); bOpy2=yy(5); kgOpy2=yy(6); qhlOpy2=yy(7);
lfOlmlf2 = alpl/(1-alpk-alpn)/thetpl/qhlOpy2/mg4py2*((1-betag*omg1)/(1-omg1) + Gamb2*mg1py2 + Gamd2*mg4py2);
Lf2 = L2*lfOlmlf2/(1+lfOlmlf2);
qhOpy2 = qhlOpy2/(L2-Lf2);

%=========
%step3
%=========
% P1 P2 Y1 Y2 Y
% ww(1) ww(2) ww(3) ww(4) ww(5)
hh=@(ww)([
    (ww(1)^(-alpk-alpg)*ww(3)^(1-alpk-alpg)/zp1/((1-Tau1)*alpk/Rk)^alpk/kgOpy1^alpg/Lf1^alpl)^(1/alpn) + ...
    (ww(2)^(-alpk-alpg)*ww(4)^(1-alpk-alpg)/zp2/((1-Tau2)*alpk/Rk)^alpk/kgOpy2^alpg/Lf2^alpl)^(1/alpn) - nbar;
%     ww(5) - ww(3)^kap*ww(4)^(1-kap);
%     kap*ww(3)*ww(1) - ww(5);
%     kap*ww(4)*ww(2) - ww(5);
    ww(5)^(1-1/epsl) - ww(3)^(1-1/epsl) - ww(4)^(1-1/epsl);
    ww(3)*ww(1)^epsl - ww(5);
    ww(4)*ww(2)^epsl - ww(5);
    varphi1*(kgOpy1)^varphig*(1-Xi1)*(1-Tau1)*alpn*(zp1*((1-Tau1)*alpk/Rk)^alpk*kgOpy1^alpg*Lf1^alpl)^(1/alpn)*qhOpy1^(btpar-1)*ww(1)^(varphig+btpar+(alpk+alpg)/alpn)*ww(3)^(varphig+btpar-(1-alpk-alpg)/alpn) - ...
    varphi2*(kgOpy2)^varphig*(1-Xi2)*(1-Tau2)*alpn*(zp2*((1-Tau2)*alpk/Rk)^alpk*kgOpy2^alpg*Lf2^alpl)^(1/alpn)*qhOpy2^(btpar-1)*ww(2)^(varphig+btpar+(alpk+alpg)/alpn)*ww(4)^(varphig+btpar-(1-alpk-alpg)/alpn)
    ]);
ww=fsolve(hh,[1,1,3,3,7],options);
P1=ww(1);P2=ww(2);Y1=ww(3);Y2=ww(4);Y=ww(5);

%=========
%step4
%=========
N1 = (P1^(-alpk-alpg)*Y1^(1-alpk-alpg)/zp1/((1-Tau1)*alpk/Rk)^alpk/kgOpy1^alpg/Lf1^alpl)^(1/alpn);
N2 = (P2^(-alpk-alpg)*Y2^(1-alpk-alpg)/zp2/((1-Tau2)*alpk/Rk)^alpk/kgOpy2^alpg/Lf2^alpl)^(1/alpn);
V1 = varphi1*(kgOpy1)^varphig*(1-Xi1)*(1-Tau1)*alpn*(zp1*((1-Tau1)*alpk/Rk)^alpk*kgOpy1^alpg*Lf1^alpl)^(1/alpn)*qhOpy1^(btpar-1)*P1^(varphig+btpar+(alpk+alpg)/alpn)*Y1^(varphig+btpar-(1-alpk-alpg)/alpn);
V2 = varphi2*(kgOpy2)^varphig*(1-Xi2)*(1-Tau2)*alpn*(zp2*((1-Tau2)*alpk/Rk)^alpk*kgOpy2^alpg*Lf2^alpl)^(1/alpn)*qhOpy2^(btpar-1)*P2^(varphig+btpar+(alpk+alpg)/alpn)*Y2^(varphig+btpar-(1-alpk-alpg)/alpn);
W1 = (1-Tau1)*alpn*P1*Y1/N1;
W2 = (1-Tau2)*alpn*P2*Y2/N2;
K1 = (1-Tau1)*alpk*P1*Y1/Rk;
K2 = (1-Tau2)*alpk*P2*Y2/Rk;
K = K1 + K2;I = delk*K;
Qf1 = (1-Tau1)*alpl*P1*Y1/Lf1;
Qf2 = (1-Tau2)*alpl*P2*Y2/Lf2;
Kg1 = kgOpy1*P1*Y1;Ig1=delg*Kg1;B1=bOpy1*P1*Y1;Qh1=qhOpy1*P1*Y1;Mug11=mg1py1/P1/Y1;Mug21=mg2py1/P1/Y1;Mug31=mg3v1/V1;Mug41=mg4py1/P1/Y1;
Kg2 = kgOpy2*P2*Y2;Ig2=delg*Kg2;B2=bOpy2*P2*Y2;Qh2=qhOpy2*P2*Y2;Mug12=mg1py2/P2/Y2;Mug22=mg2py2/P2/Y2;Mug32=mg3v2/V2;Mug42=mg4py2/P2/Y2;
C1 = btpar*(1-Xi1)*W1;
C2 = btpar*(1-Xi2)*W2;
Lh1 = (1-btpar)*(1-Xi1)*W1/Qh1;
Lh2 = (1-btpar)*(1-Xi2)*W2/Qh2;
Ls1 = L1 - Lf1 - N1*Lh1;
Ls2 = L2 - Lf2 - N2*Lh2;


%% Here we define the steady state values of the endogenous variables of the model.
ys=[C1
Y1
Ig1
Kg1
K1
B1
V1
N1
P1
W1
Qh1
Qf1
Lh1
Lf1
Ls1
Mug11
Mug21
Mug31
Mug41
C2
Y2
Ig2
Kg2
K2
B2
V2
N2
P2
W2
Qh2
Qf2
Lh2
Lf2
Ls2
Mug12
Mug22
Mug32
Mug42
Y
K
I
Rb
Rk
Gamb1
Gamd1
Xi1
Tau1
L1
Z1
Gamb2
Gamd2
Xi2
Tau2
L2
Z2];
