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

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;
%%
A = 1;
Pi = Pi_bar;
Rb = 1/betta;
L = s;
N = 1;
Nw = wn*N;
Nf = (1-wn)*N;
U = s;
MatC = chi*A*(L/U)^iot;
J_by_JJ   = ( (1 - thetta*Pi^(epsilon_p-1))/(1-thetta) )^(1/(1-epsilon_p));
MC = J_by_JJ * ((epsilon_p-1)/epsilon_p) * ( 1 - thetta*betta * Pi^epsilon_p ) / ( 1 - thetta * betta * Pi^(epsilon_p-1) ); 
Disp = ( (1-thetta)*J_by_JJ^(-epsilon_p) ) / ( 1- thetta * Pi^epsilon_p );
options = optimset('TolFun',1e-9,'TolX',1e-9,'MaxIter', 5000, 'MaxFunEvals', 5000, 'Display','on');
H =  fsolve(@(x) b + (x^(1+varphi_h)/(1+varphi_h))+ xi_w/(1- xi_w)*MatC - ((1-alp)*MC*A*(x*Nw)^(1-alp))/Nw/ww - MatC/ww*(betta_e*(1-s)-1)/(( 1 + (1 - betta_e*xi_b*Rb)*xi_e/Rb)*(cy*A * (Nw*x)^(1-alp)/Disp)^(-sigma_c)),0.7,options);
W = (b + (H^(1+varphi_h)/(1+varphi_h))+ xi_w/(1- xi_w)*MatC)/H;
Wc = ww*W;
Wg = (1-ww)*W;
Yw = A * (Nw*H)^(1-alp);
Y = Yw/Disp;
G_bar = gy*Y;
G = G_bar;
C_e = cy*Y;
lambda_e = C_e^(-sigma_c);
lambda_c = lambda_e - betta_e*xi_b*lambda_e(+1)*Rb;
C = Y - C_e - G - MatC*L;
Co = wc*C;
Cp = (1-wc)*C;
lambda_b = C^(-sigma_c);
J = ( epsilon_p/(epsilon_p-1) ) * Yw * MC/ (1- thetta * betta * Pi^epsilon_p); 
JJ = Yw/ (1- thetta * betta * Pi^(epsilon_p-1));
YY=1;
CC=1;
NN = 1;
UU = 1;
HH = 1;
RR=1;
PiPi=1;
dyobs= 0;
dcobs= 0;
dnobs= 0;
duobs= 0;
dhobs= 0;
pinfobs= 0; 
shadow_robs= 0;

%% end own model equations

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

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