function [ys,params,check] = dsge1_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = dsge1_steadystate(ys,exo,M_,options_)
% computes the steady state for the dsge1.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
%   - 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] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright (C) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

%% 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
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

%% Step 2: Enter parameter restrictions
% if ca2*cb2<1 % parameter viloates restriction
%     check=1; % set failure indicator
%     return; % return without updating steady states
% end

%% Step 3:  Enter model equations here

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

% the steady state computation follows FVRR (2006), section 4.1
 c = ((c_t)^ca1 * (c_n)^(1-ca1))/(ca1^ca1*(1-ca1)^(1-ca1));
1-ca1/((p_n)*(c_n)= ca1/(e)*(c_t);
    cb2/((p)*(c)) = cb1/(m)  + cb3 * cb2/(p(+1)*c(+1));
         ca2/cb2 = w*(1-lab)/(p)*(c);
              p_t = e;
                p = (p_t)^(ca1)*(p_n)^(1-ca1);
              c_t = (x)*(k_t)^(cf1)*(lab_t)^(1-cf1) - g_t + (t-u)*(q);
              c_n = (z)*(k_n)^(cf2)*(lab_n)^(1-cf2) - (g_n);
                r = cf1*(p_t)*(x)*(k_t)^(cf1-1) *(lab_t)^(1-cf1);
                r = cf2*(p_n)*(z)*(k_n)^(cf2-1) *(lab_n)^(1-cf2);
                w = (1-cf1)*(x) * p_t *(k_t)^(cf1) *(lab_t)^(-cf1);
                w = (1-cf2)*(z)* p_n *(k_n)^(cf2) *(lab_n)^(-cf2);
                k = k_t + k_n;
              lab = lab_t + lab_n;
                y = y_t + y_n;
              y_t = (x)*(k_t)^(cf1) *(lab_t)^(1-cf1);
              y_n = (p_n/e) * (z)*(k_n)^(cf2) *(lab_n)^(1-cf2); 
                g = ch2*g_n + (1-ch2)*g_t;
              g_n = ch2*g;
              g_t = (1-ch2)*g;
              

%% 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 %auxiliaryvariables are set automatically
%      eval([ 'ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
%  end
     
%% Step 4: Update parameters and variables

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

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