function [ys,check] = li_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

beta1=1.5;
h=0.5;
gamma_l=2;
kappa=4;
theta_w=0.6;
chi_w=0.5;
alpha=0.45;
theta_p=0.6;
chi=0.5;
gamma_y=2;
gamma_pi=2;
Lambda_mu1=0.62;
yta=10;
gamma_w=0.5;
rho_d=0.5;
rho_phi=0.5;

beta = 1.0/(beta1/100.0+1.0);
Lambda_mu=Lambda_mu1 / 100;
muss=exp(Lambda_mu);
piss=1.072;
lss=0.25;
zss=exp(0.0238);

%gamma_1=r;
delta = 1.0-(1.0-0.0678)*zss*muss; %depreciation rate
r = zss*muss/beta-(1.0-delta); %gamma1
R = piss*zss/beta; %s.s. nominal interest rate
Ass = exp(log(zss)*(1-alpha)-log(muss)*alpha); %s.s. growth rate of neutral technology shock

pi_s = ((1.0-theta_p*piss^((1.0-chi)*(eps-1.0)))/(1.0-theta_p))^(1.0/(1.0-eps)); %s.s. ratio of optimal price to aggregate price
Delta_p = (1.0-theta_p)*pi_s^(-eps)/(1.0-theta_p*piss^((1.0-chi)*eps)); %s.s. relative price dispersion
bigpsi = ((eps-1.0)/eps)*pi_s*(1.0-beta*theta_p*piss^((1.0-chi)*eps))/(1.0-beta*theta_p*piss^((1.0-chi)*(eps-1.0))); %s.s. real marginal cost
w = (1.0-alpha)*(bigpsi*(alpha/r)^alpha)^(1.0/(1.0-alpha)); %s.s. real aggregate wage
k = ld*(alpha/(1.0-alpha))*(w/r)*zss*muss; %s.s. capital-labor ratio
F = Delta_p/5.873-(Ass/zss)*(k/ld)^(alpha-1.0);
eps=x;
eps0=1.5;
[x,Fval,exitflag] = fsolve(@(x) F,eps0,optimset('Display','off'),);

pi_w = ((1.0-theta_w*piss^((1.0-chi_w)*(yta-1.0))*zss^(yta-1.0))/(1.0-theta_w))^(1.0/(1.0-yta)); %s.s. ratio of optimal wage to aggregate wage
Delta_w = (1.0-theta_w)*pi_w^(-yta)/(1.0-theta_w*zss^yta*piss^((1.0-chi_w)*yta)); %s.s. relative wage dispersion
ld = lss/Delta_w; %s.s. effective labor service
w_s = w*pi_w; %s.s. real optimal wage
%//Omega = k_Ld; %auxiliary variable for s.s. capital-labor ratio
%//k_tilde = Omega*Ld; %s.s. capital stock
%//i_tilde = i_k_tilde*k_tilde; %s.s. investment
%//y_d_tilde = k_tilde/k_y_tilde; %s.s. output
c=yd-i;
%//c_tilde = y_d_tilde*(1.0-i_tilde/y_d_tilde); %s.s. consumption
lam = (1.0-h/zss)^(-1.0)*c^(-1.0); %s.s. shadow price
nu_m = lam*(1.0-beta/(zss*piss)); %s.s. ratio related to money demand 
%//LHS = (1.0-beta*theta_w*zss^(yta*(1.0+gamma_l))*piss^(yta*(1.0-chi_w)*(1.0+gamma_l)))/(1.0-beta*theta_w*zss^(yta-1.0)*piss^((1.0-chi_w)*(yta-1.0)));
%//psi = (LHS*((yta-1.0)/yta)*wss_s*lamss)/(piss_w^(-yta*gamma_l)*ldss^gamma_l); %weight on labor service
i=k*0.0678;
k=y*5.873;
omega=zss*piss;

rss=r;
Rss=R;
piss_s=pi_s;
Deltass_p=Delta_p;
bigpsiss=bigpsi;
wss=w;
kss=k
piss_w=pi_w;
Deltass_w=Delta_w;
ldss=ld;
wss_s=w_s;
css=c;
lamss=lam;
nu_mss=nu_m;
omegass=zss*piss;


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