function [ys,check]=Benin_steadystate(ys, exe);

global M_ lgy_ 
global theta1 theta2 theta3 bet deltaP deltaG rhoZ rhoGC rhoGI khi phiGC phiGI sigmaZ sigmaGC sigmaGI
   
if isfield(M_ lgy_.'param_nbr') == 1
NumberOfParameters = M_ lgy_.param_nbr;
for i = 1:NumberOfParameters
  paramname = deblank(M_ lgy_.param_names(i,:));
  eval([ paramname ' = M_ lgy_.params(' int2str(i) ');']);
end
check = 0;
end

% Model equations to be entered here

global theta1 theta2 theta3 bet deltaP deltaG rhoZ rhoGC rhoGI khi phiGC phiGI sigmaZ sigmaGC sigmaGI

% theta1=0.202;   Private capital share
% 
% theta2=0.087;   Public capital share
% 
% theta3=0.711;   labor (wage) share
%
% bet=0.9506;     discount factor
%
% deltaP=0.097;   depreciation rate of private  capital
%
% deltaG=0.054;   depreciation rate of public capital
%
% rhoZ=0.832;     Productivity shock persitance (AR(1) parameter of TFP)
%
% rhoGC=0.945972; government consumption spending persistance parameter (AR(1) term)
%
% rhoGI=0.630941; government investment spending persistance parameter (AR(1) term)
%
% sigmaZ=0.033;   standard deviation of technical shock
%
% sigmaGC=0.080;  standard deviation of Gvt consumption spending shock
%
% sigmaGI=0.411;  standard deviation of Gvt nvestment spending shock
%
% khi=0.8;        Utility weight of leisure
%
% phiGC=0.018833; Correlation between GDP and gvt consumption expenditure
%
% phiGI=0.017672; Correlation between GDP and gvt consumption expenditure

x0 = [0.149; 0.1401; 0.0737; 0.0776; 0.374; 0.0999; 0.0097; 0.4757; 0.0035;  0.0257; -0.0292; 1; 0: 0; 0];
[x,fval,exitflag] = fsolve(@Benin_steadystate, x0); % Solving the nonlinear system
if exitflag < 1
    check=1;
end

y = x(1);
cp = x(2);
cg = x(3);
kp = x(4);
kg = x(5);   
invp = x(6);
invg = x(7);
h = x(8);
w  = x(9);
r = x(10);
pb= x(11);
z = x(12);
tau_C=x(13);
tau_H=x(14);
tau_K=x(15);

z=1;
r=(1-bet+bet*deltaP)/(bet*(1+tau_K));
h=(((1+tau_C)/(1-tau_H)*theta3)*(1-(phiGI/(1-rhoGI)+phiGC/(1-rhoGC)+theta1*deltaP/r)))^(1/1+khi);
y=(h^theta3*(theta1/r)^theta1*(phiGI/(1-rhoGI)*deltaG)^theta2)^(1/(1-theta1-theta2));
cp=((1-tau_H)*theta3*y)/((1+tau_C)*(h)^(1-khi));
w=theta3*y/h;
kp=theta1*y/r;
invp=theta1*y*deltaP/r;
kg=(phiGC*y)/((1-rhoGC)*deltaG);
invg=(phiGI*y)/(1-rhoGI);
cg=(phiGC*y)/(1-rhoGC);
pb=tau_C*cp+(theta1*tau_K+theta3*tau_H)*y-cg-invg;
tau_H=0;
tau_C=0;
tau_K=0;

for iter = 1:length(M_ lgy_.params)
  eval([ 'M_ lgy_.params(' num2str(iter) ') = ' M_lgy_.param_names(iter,:) ';' ])
end

if isfield(M_ lgy_,'param_nbr') == 1


NumberOfEndogenousVariables = M_ lgy_.endo_nbr; %auxiliary variables are set automatically
ys=zeros(NumberOfEndogenousVariables,1);
for ii = 1:M_ lgy_.orig_endo_nbr
  varname = deblank(M_ lgy_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end

end