function [ys,check] = dynamics_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval(['global ' paramname]); % The line is added by YB. Note global variables are expected to be claimed before given values.
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%% Enter model equations here
%% Solve the steady state
q=1-1/psiG;
A=1;    % A is dropped in the following calculation
qc=1;    % qc is dropped in the following calculation
fn_w=@(w) fn_main3(w,q);
w=fsolve(fn_w,1.5);
[yy,themat,Vaug,k1,alpmat,epsmat,h,iter]=fn_main3(w,q);

pi=1;
chi=0;
R=1/rhoG;
i=0;
zet=( bet*q/w )^( 1/(1-bet) );
y=0;
k=0;
for j=1:J
    i=i+themat(j)*alpmat(j)*k1*( 1-(1-del)^j );
    y=y+themat(j)*zet^bet*k1^( alp/(1-bet) )*(1-del)^( alp/(1-bet)*(j-1) );
    k=k+themat(j)*k1*(1-del)^(j-1);
end
c=y-i;
lamd=(1-b*rhoG)*( (1-b)*c-h0*h^(1+h1)/(1+h1) )^(-sig);

%% display the results
disp('The no. of the iteration is:')
iter
disp('check whether one equation is equal to 0:')
lamd*w-h0*h^h1*( (1-b)*c-h0*h^(1+h1)/(1+h1) )^(-sig)

% ajsustment costs relative output
adjcost=0;
for j=1:J
    adjcost=adjcost+themat(j)*w*fn_dist_cost(alpmat(j));
end
disp('The ratio of adjustment cost to investment is (%):')
adjcost/i*100

for j=1:J
    eval(sprintf('the_%d = themat(j);',j))
end
for j=1:(J-1)
    eval(sprintf('alp_%d = alpmat(j);',j))
	eval(sprintf('V_%d = Vaug(j+1);',j))
	eval(sprintf('xi_%d = epsmat(j);',j))
end
eval(sprintf('V_%d = Vaug(1);',0))

eta=(1-bet)*bet^( bet/(1-bet) )*(q*A)^( 1/(1-bet) )*w^( -bet/(1-bet) );
R_ss=R;
y_ss=y;
kmat=zeros(J,1);
for j=1:J
    kmat(j)=(1-del)^(j-1)*k1;
	eval(sprintf('k_%d = kmat(j);',j))
end

%% end own model equations

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

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