function [ys,check] = try4_steadystate(ys,exo)
% function [ys,check] = try4_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;


%% Enter model equations here

options=optimset(); % set options for numerical solver

% the steady state computation follows FV (2006), section 4.1
S=S_bar;
P_f=P/S;
I=1/bet;
I_f=I;

function F = myfun(x)
F = [x(1)^(1-theta)*tau^(theta-1)*(k/(k-theta+1))^2+(x(1)^(-k)*(k/(k-theta+1))^(k/(theta-1))*((theta-1)/(k-theta+1))+(bet-bet*delta-1)*fE/fX);
     x(2)-((S/P)^(theta-1)*((theta*k*fX*I)/(k-theta+1)-fE*((1-bet)/bet)/((tau*zd/x(1))^(theta-1)+(1/x(1))^k*(k/(k-theta+1))^(k/(theta-1))))^(1/(1-theta)))];
[x,fval] = fsolve(@myfun,x0,options)
 zx=x(1);
 PX=x(2);

ND=((S*PX/P)^(theta-1))/((tau*zd/zx)^(theta-1)+(1/zx)^k*(k/(k-theta+1))^(k/(theta-1)));
PD=S*zx*PX/(tau*zd);
W=S*PX*(theta-1)*zx/(theta*tau);
C=W*(L+ND*fE*(1-bet)/bet)/P;
NE=delta*ND;
NX=ND*(zx)^(-k)*(k/(k-theta+1))^(k/(theta-1));
dd=(PD/P)^(1-theta)*C*P/theta;
dx=(S*PX/P)^(1-theta)*C*P/theta-W*fX;
v=W*fE;
d=dd+(NX/ND)*dx;
zx_f=zx;
PX_f=PX;
ND_f=ND;
PD_f=PD;
W_f=W;
C_f=C;
NE_f=NE;
NX_f=NX;
dd_f=dd;
dx_f=dx;
v_f=v;
d_f=d;





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