function [ ys,check ] = Neoclassical_Model_SP_Hidden_Chi_NEW_steadystate( ys,exo )
% function [ ys,check ] =Neoclassical_Model_SP_Hidden_Chi_0_5_steadystate( ys,exo )
% computes the steady state for the Neoclassical_Model_SP_Hidden_Chi_0_5.mod and uses a numerical
% solver to do so
% Inputs:
% - ys [vector] vector of initial values for the steady state of the
% endogenous vaiables
% - exo [vector] vector of values for the exogenous variables
%
% Output:
% - ys  [vector] vector of steady state values for the endogenous variables
% - check  [scalar] set to 0 if steady state computation worked and to 1 if
% not (allows to impose 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;


%% Enter model equations here
%Alternative Steady State solver through numerical methods
%syms c k l

%F1 = psii*c*(l^chii) - (1-alppha)*(k^alppha)*(l^(-alppha));
%F2 = 1/betta - (1 + alppha*(k^(alppha-1))*l^(1-alppha) - delta);
%F3 = c - (k^alppha)*(l^(1-alppha)) + delta*k;

%X = vpasolve(F1, F2, F3);

%c = double(X.c);
%k = double(X.k);
%l = double(X.l);

% Set of Steady State Equations (pinned down by hand - deterministic SS)
% k/l = (((1/betta - 1 + delta)/alppha)^(1/(alppha -1)))
% l = [((1-alppha/psii)*(k/l)^alppha)/((k/l)^alppha - delta*(k/l))]^(1/1+chii);
l = ((1-alppha/psii)*((((1/betta - (1 - delta))/alppha)^(alppha/alppha-1)))/(((1/betta - (1 - delta))/alppha)^(alppha/alppha-1))-delta*(((1/betta - (1 - delta))/alppha)^(1/alppha-1)))^(1/1+chii);
% Alternative (previous) way to pin down l:
%l = (((1 - alppha)/psii)*(((1/betta) - 1 + delta)/((1/betta) - 1 + delta - alppha*delta)))^(1/1+chii);
k = (((1/betta - 1 + delta)/alppha)^(1/(alppha -1)))*l;
y = (k^alppha)*(l^(1-alppha));
invest = delta*k;
c = y - invest;
y_l = y/l;
z = 0;
log_y=log(y);
log_invest=log(invest);
log_c=log(c);
log_l=log(l);
log_k=log(k);
log_y_l=log(y_l);

%% 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


end