function [ys,check] = drh1nk_steadystate(ys,exo)
% function [ys,check] = alpha_steadystate(ys,exo)
% computes the steady state for the alpha.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;

%% 
alpha = 1/3;
betaa = 0.9;
delta = 0.025;
psi = 0.1;
theta = 1.1;
epsilon = 0.5;
pii_ss = 0.01;
tau = 0.1;
rho = 0.9;


z = 1;
p2 = 1;
p1 = 1;
mc = 1/theta;
r = 1/betaa-1+delta;
R = pii_ss/betaa;
pii = pii_ss

%%
x0 = [1,0.9];
[xx,fval]=fsolve(@nkss,x0,optimoptions('fsolve','TolFun',1e-5, 'MaxIter', 1e15),alpha,betaa,delta,psi,theta,epsilon,pii_ss,tau);
k = xx(1);
c = xx(2);

lambda = 1/c;
i = delta*k;
w = (mc/(1/alpha)^alpha/(1/(1-alpha))^(1-alpha)/r^alpha)^(1/(1-alpha));
h = (lambda*w/tau)^(1/psi);
y = p1^(theta/(theta-1))*k^alpha*h^(1-alpha);
fd = lambda*y/(1-betaa*epsilon);
kd = lambda*y*theta*mc/(1-betaa*epsilon);

%% 

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
