function [ys,params,check] = example_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;

%% Enter model equations here
load('stst_guess.mat');
run solve_SS

fL_2=theta_L^(1-eta)*muL_2;
fH_2=theta_H^(1-eta)*muH_2;
fL_1=theta_L^(1-eta)*muL_1;
fH_1=theta_H^(1-eta)*muH_1;
f_L=fL_2+fL_1;
f_H=fH_2+fH_1;
q_L=f_L/theta_L;
q_H=f_H/theta_H;

uuuH=d*ff_H/(d+f_H);
eeeH_1=(fH_1*uuuH)/d;
eeeH_2=ff_H-uuuH-eeeH_1;

uuuL=d*ff_L/(d+f_L);
eeeL_1=(fL_1*uuuL)/d;
eeeL_2=ff_L-uuuL-eeeL_1;

v_L=theta_L*uuuL;
v_H=theta_H*uuuH;
M_L=((uuuL)^eta)*(v_L^(1-eta));
M_H=((uuuH)^eta)*(v_H^(1-eta));

Y_2=(alpha_2*((eeeL_2)^((gggamma-1)/gggamma))+(1-alpha_2)*((eeeH_2)^((gggamma-1)/gggamma)))^(gggamma/(gggamma-1));
Y_1=(alpha_1*((eeeL_1)^((gggamma-1)/gggamma))+(1-alpha_1)*((eeeH_1)^((gggamma-1)/gggamma)))^(gggamma/(gggamma-1));
P= (p1*(p1/a)^(-sssigma)+p2*(p2/(1-a))^(-sssigma))^(1/(1-sssigma));

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
