function [ys,check] = modelv01_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
check=0;
%% Enter model equations here
cgamma=-cgammaM;
UpsilonP=1;
PinBar=1;
PiC=PiC_SS;
Pnbar=1;
Pobar=1;
Pin=PiC_SS;
Pio=PiC_SS;
M=1;
MC=(cepsilon-1)/cepsilon;
R=PiC_SS/cbeta;
OYratio=cphio*MC;
L=L_SS;
Y=(L^cphin*OYratio^cphio)^(1/(1-cphio));
O=OYratio*Y;
C=Y-O;
Cn=(1-cpsi)*C;
Co=cpsi*C;
F1=MC*Y/(1-cbeta*cksin);
F2=Y/(1-cbeta*cksin);
W=cphin*MC*Y/L;
cchi=W/((LMax-L)^(-csigmaL)*C^csigmaC);
V=(C^(1-csigmaC)/(1-csigmaC)-cchi*(LMax-L)^(1-csigmaL)/(1-csigmaL))/(1-cbeta);
VK=V;
VExp=1;
V_SS=V;
Y_SS=Y;
logR=log(R);
logY=0;
logPiC=0;
DIV=Y-W*L/Pnbar-Pobar*O/Pnbar;
QE=cbeta*DIV/(1-cbeta);
%% 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
