function [ys,params,check] = tank3_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] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% 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

% initialize indicator
check = 0;

options=optimset('display','off'); % set options for numerical solver
options = optimoptions('fsolve','Display','iter');

rr=1/bbeta;
pi=1;
r=rr*pi;
er=1;
eg=1;
mc=(eeta-1)/eeta;
y=1;
c=y;
n=1/3;
a=y/n;
w=a*mc;
c_r=c;
c_nr=c;
n_r=n;
n_nr=n;
g_nr=c_nr-w*n_nr;
g=(1-kka)*g_nr;
tax=0.282;
b=(tax-g)/(rr-1);
x_1=y*mc/(1-zzeta*bbeta);
x_2=y/(1-zzeta*bbeta);
p_reset=1;
llambda=c_r^(-ssigma);



kka_n=w/(n_r^(ggamma)*c_r^(-ssigma));


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