function [ys,params,check] = RBC_credit_det_steadystate(ys,exo,M_,options_)
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] 0 if steady state computation worked and to
%                        1 of not (allows to impose restrictions on parameters)

%% Step 0: initialize indicator and set options for numerical solver
check = 0;
options = optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

%% Step 1: read out parameters to access them with their name
for ii = 1:M_.param_nbr
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
end
params=M_.params;

%% Step 2: Check parameter restrictions
%if etaC*etaL<1 % parameter violates restriction (here it is artifical)
%    check=1; %set failure indicator
%    return;  %return without updating steady states
%end

%% Step 3: Enter model equations here

% What does setting to 1 would mean?
%z = epsz/(1-RHOZ);
z = 1;

A = 1;

p = 1;

m = 1;

% use different values 
g = 1;

x = 1;

if z <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% r = 0.035
r = 1/BETA + DELTAK -1;

if r <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% K_NH 740.35 re-check the -values 
K_NH = (r/(ALPHA*exp(z)))^(1/ALPHA-1);

if K_NH <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% Y_NH = 70.69
Y_NH = exp(z)*K_NH^ALPHA;

if Y_NH <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% IV_NH = 1.767
IV_NH = DELTAK*K_NH;

if IV_NH <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% w = 45.95
w = (1-ALPHA)*exp(z)*K_NH^ALPHA;

if w <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% e = 0.0005 (when omega = 10, e = 0.069)
e = (w/(OMEGA*PHIE))^(1/(PHIE - 1));

if e <= 0.05
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% n = 0.026 use fsolve here n0 = 0.35
n = (1/w)*(1/BETA - (1-DELTAH) - OMEGA*e^(PHIE));

if n <= 0.1
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

% might need to use fsolve asthe value is negative
h = x/(DELTAH - OMEGA*e^(PHIE));

%X_NH = 1/(n*(DELTAH - OMEGA*e^(PHIE)));

if h <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

X_NH = x/(n*h);

C_NH = Y_NH - X_NH - IV_NH;

if C_NH <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end


l = 1 - n - e;

if l <= 0
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

c = ((BETA*w*h*l^(GAMMA))/g*A)^(1/SIGMA);

if c/n*h == C_NH
    check = 1; % set failure indicator
    return;    % return without updating steady states
end

p_hat = 1/c;

M = p/p_hat;

m_hat = m/M;

y  = Y_NH*n*h;
iv = IV_NH*n*h;
k  = K_NH*n*h;



% The labor level
%if etaC == 1 && etaL == 1
%    % Closed-form solution for labor
%    L = gam/pssi*C_L^(-1)*W/(1+gam/pssi*C_L^(-1)*W);
%else
%     No closed-form solution use a fixed-point algorithm
%    L0 = 1/3;
%    [L,~,exitflag] = fsolve(@(L) pssi*(1-L)^(-etaL)*L^etaC - gam*C_L^(-etaC)*W, L0,options);
%    if exitflag <= 0
%        check = 1; % set failure indicator
%        return     % return without updating steady states
%    end
%end



%% Step 4: Update parameters and variables
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

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