function [ys,check] = model_steadystate(ys,exe)
global M_ lgy_

if isfield(M_,'param_nbr') == 1
NumberOfParameters = M_.param_nbr;
for i = 1:NumberOfParameters
paramname = deblank(M_.param_names(i,:));
eval([ paramname ' = M_.params(' int2str(i) ');']);
end
check = 0;
end

%%%% Model equations to be entered here


as=1;
rs = ((1+thet)/bet) -1 + delt - (thet/PIS);
ws = (1-alph)*((rs/alph)^(alph/(alph-1)));
IS = PIS/bet;
f = @(l2) (1+et*(PIS-1))*(((PIS*et - et*bet + bet)*ch/(bet*ws))^(-1/sigm))*(l2^(-gamm/sigm)) - ws*l2;
l2s = fzero(f,1);
c2s = (((PIS*et - et*bet + bet)*ch/(bet*ws))^(-1/sigm))*(l2s^(-gamm/sigm));
m2s = et*PIS*c2s;
g = @(l1) ws*l1 + (rs-delt)*(l1+l2s)*((rs/alph)^(1/(alph-1)))*(rs-delt+(thet/PIS)*(1-IS)) + ((PIS-IS)/PIS)*m2s - (ws*(l1^(-gamm))/ch)^(1/sigm);
l1s = fzero(g,1);
c1s =  (ws*(l1s^(-gamm))/ch)^(1/sigm);
ks = (l1s+l2s)*((rs/alph)^(1/(alph-1)));
m1s = thet*ks;
bs = ((PIS-1)/(IS-PIS))*IS*(m1s+m2s);
ls = l1s + l2s;
lambds = ((ws/(1-alph))^(1-alph))*(rs/alph)^alph;
ys = (ks^(alph))*(ls^(1-alph));
c1 = 0;
c2 = 0;
l1 = 0;
l2 = 0;
m1 = 0;
m2 = 0;
k = 0;
b = 0;
r = 0;
w = 0;
ii = 0;
lambd = 0;
l = 0;
pii = 0;
a = 0;
ea = 0;
epi = 0;
nu = 0;
y = 0;

for iter = 1:length(M_.params)
    eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

if isfield(M_,'param_nbr') == 1

if isfield(M_,'orig_endo_nbr') == 1
NumberOfEndogenousVariables = M_.orig_endo_nbr;
else
NumberOfEndogenousVariables = M_.endo_nbr;
end
ys = zeros(NumberOfEndogenousVariables,1);
for i = 1:NumberOfEndogenousVariables
varname = deblank(M_.endo_names(i,:));
eval(['ys(' int2str(i) ') = ' varname ';']);
end
else
ys=zeros(length(lgy_),1);
for i = 1:length(lgy_)
ys(i) = eval(lgy_(i,:));
end
check = 0;
end
end