function [ys,check] = RBC_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([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%steady states
q_ss=1; 
A_ss=1;
R_ss=1/betta; 
rk_ss=1/betta+delta-1; 
K_ss=(rk_ss/alppha)^(1/(alppha-1))*N_ss;
Y_ss=K_ss^alppha*N_ss^(1-alppha);
w_ss=(1-alppha)*Y_ss/N_ss; 
I_ss=delta*K_ss; 
C_ss=Y_ss-I_ss; 
lambda_ss=(C_ss*(1-hab))^(-siggma); 
psi=w_ss*lambda_ss*N_ss^(-eta);

prod_ss=Y_ss/N_ss; 

lambda_t=log(lambda_ss); 
C_t=log(C_ss); 
R_t=log(R_ss); 
rk_t=rk_ss; 
w_t=log(w_ss); 
N_t=log(N_ss); 
Y_t=log(Y_ss); 
A_t=log(A_ss); 
K_t=log(K_ss); 
I_t=log(I_ss); 
SDF_t=betta; 
q_t=log(q_ss);
prod_t=log(prod_ss);
AA_t=0;

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