function [ys,params,check] = Test_steadystate(ys,exo,M_,options_)

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

%%

% global betta alfa epsilon delta ies chi csi iota theta bP_Y piss ps deltax D0 D1 D2  pl G ...
% levss nuss spss zeta levGss kGss kss kappaI phipi phiy rhor rhoc rhoa rho_xi  erow rr ...
% kappaP rss bPss levBss bstar  spss spGss spBss omega sig alpha0 muGss muss kappaFG psi  theta1 theta2 taubar mueps mugeps

load par_ce.mat;

tau = exo(strmatch('tau',M_.exo_names,'exact'));
mu_eps = exo(strmatch('mu_eps',M_.exo_names,'exact'));
muG_eps = exo(strmatch('muG_eps',M_.exo_names,'exact'));

taubar = tau;
mueps = mu_eps;
mugeps = muG_eps;


load x0                      % Initial guess
spG=spss;                                         % green spread            
spB=spss;                                         % brown spread
rG=rr+spG;                                      % loan rate for green bonds
rB=rr+spB;                                      % loan rate for brown bonds
rkG=rG-(1-delta);                               % rental rate of green capital
rkB=rB-(1-delta);                               % rental rate of brown capital
pI=(epsilon-1)/epsilon; 

%taubar=0.1;
options = optimoptions('fsolve','MaxFunEvals',300000,'MaxIter',30000,'TolFun',1e-15);
x = fsolve(@(xx) FindSS(xx,alfa,omega,delta,csi,rkB,rkG,pI,G,.2,pl,taubar,ps,zeta,deltax,erow,D0,D1,D2),x0,options);

y=x(1);
pB=x(2);
k=x(3);
e=x(4);

g=G*y;
kG=.2*k;
kB=k-kG;
    yB = (rkB * kB) / (alfa * pB);
    for i = 1:1000  % Iterative approximation (Newton’s method could be used here)
        yB_new = (rkB * kB) / (alfa * (pB - taubar *(1-ps)* (yB^(1-ps))));
        if abs(yB_new - yB) < 1e-20  % Convergence condition
            break;
        end
        yB = yB_new;
    end
X=(e+erow)/deltax;
A=(1-D0-D1*X-D2*X^2);
pG=(1/(1-zeta)*(pI^(1-csi)-zeta*((pB)^(1-csi))))^(1/(1-csi));
yG=(1-zeta)*(pG/pI)^(-csi)*y;   
Lb=(yB/(A*kB^(alfa)))^(1/(1-alfa));
Lg=(yG/(A*kG^(alfa)))^(1/(1-alfa));
iB=delta*kB;
iG=delta*kG;
i=iG+iB;  
c=y-i-g; 
wG=(1-alfa)*pG*yG/Lg;
L=(Lb^(1+pl)+Lg^(1+pl))^(1/(1+pl));

wB=(1-omega)/omega*c/((1-L)*(Lb^(1+pl)+Lg^(1+pl))^(pl/(1+pl))*Lb^(-pl));
          
%% Other steady-state values
lev=levss;
         
b=kB+kG;                     % total bonds
bG=kG;                       % green bonds
bB=kB;                       % brown bonds
bP=bP_Y*4*y;                 % bonds held by the public sector
muG=muGss + mugeps;
mu=muss+mueps;
%mu=bP/b;
bPG=mu*muG*b;                  % green bonds held by the public sector
bPB=bP-bPG;                  % bonds held by the public sector
bF=b-bP;                     % bonds held by the financial sector
n=bF/lev;                     % net worth 
bFG=bG-bPG;                  % green bonds held by the financial sector
bFB=bB-bPB;                  % brown bonds held by the financial sector
levG=(kG-muG*bP)/(bF/lev);
%levB=lev-levG;

%% end own model equations

a=1; 
M=betta;
nu=nuss;
muG=muGss + mugeps;
r=piss/betta;
r_temp=piss/betta;
V = ((c^omega * (1-L)^(1-omega))^(1-1/ies) -(alpha0*X)^(1-1/ies))^(1/(1-1/ies)) ;
CE = V^(1-sig); 
%D=n-bFG+bFB;  
qB=1;
qG=1;
pi=piss;
Xib=1;
Xig=1;
rBE=rB;
rGE=rG;
spGB=0;
sg=bG/b;
skG=kG/k;
y_gap=0;
dr=0;
L_CB=piss^2;


params=NaN(NumberOfParameters,1);
for iter = 1:NumberOfParameters
    param_value = eval(M_.param_names{iter});
    params(iter) = param_value;
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
