function [ys,params,check] = SOE_one_sector_steadystate(ys,exe,M_,options_)

%% Load model parameters

NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = M_.param_names{i,:};                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end         
check = 0;
%%

varrho_0=0.7;
logit_varrho_0=log(varrho_0/(1-varrho_0));
options = optimset('TolFun',1e-9,'TolX',1e-9,'MaxIter', 5000, 'MaxFunEvals', 5000, 'Display','on');
[logit_varrho] = fsolve(@(logit_varrho)ss_fun(logit_varrho,M_),logit_varrho_0,options);
varrho=exp(logit_varrho)/(1+exp(logit_varrho))
Pi = Pi_bar;
PiF = Pi_bar;
Pi_star = Pi_star_bar;
PiS = PiF/Pi_star;
P_star=P_star_bar;
phi = ((Pi/betta) / R_star_bar) / PiS;
J_by_JJ   = ( (1 - xi*Pi^(zzeta-1))/(1-xi) )^(1/(1-zzeta));
MC = J_by_JJ   * ((zzeta-1)/zzeta) * ( 1 - xi*betta * Pi^zzeta ) / ( 1 - xi * betta * Pi^(zzeta-1) ); 
Disp = ( (1-xi) * J_by_JJ  ^(-zzeta) ) / ( 1- xi * Pi^zzeta );
h_by_K = ( ( 1/betta - (1-delta) ) / ( MC*(1-alp) ) )^(1/alp);
WP = alp * MC * h_by_K^(alp-1);
Yw =  H_bar * h_by_K^(alp-1);
Y = Yw/Disp;
BF =-Y * log(phi) / phiB;
TB = BF * ( 1/(phi*R_star_bar) - PiS/Pi );
C = Y+YO -  delta*H_bar/h_by_K -  gy*Y - TB;
EX = Y - delta*wI*H_bar/h_by_K - C*wC - gy*Y;

I_star_bar = isexp*EX;
C_star_bar = csexp*EX;
Gov_bar = gy*Y;
HRoT=(1-varrho)/(1-varrho*chi);
CRoT=WP*HRoT;
A = 1;
MS = 1;
PH_by_P = 1;
PF_by_P = 1;
PH_by_PI = 1;
PF_by_PI = 1;
PI_by_P = 1;
SX = 0;
SXp = 0; 
Lambda = betta;
Terms = 1;
PiH = Pi_bar;
Q = 1;
R = Pi/betta;
%CC = HC*WP;
K = Yw/( h_by_K^alp );
H = K*h_by_K;
HR = H/(1-lambda) - HRoT*lambda/(1-lambda);
I = delta*K;
IH = wI*I;
IF = (1-wI)*I;
Gov = gy*Y;
CH = C*wC;
CF = (1-wC)*C;
CR = C/(1-lambda) - CRoT*lambda/(1-lambda);
cr = log( CR );
J = ( (zzeta/(zzeta-1))*Y*MC ) / (1 - xi * betta * Pi^zzeta ); 
JJ = Y / ( 1 - xi * betta * Pi^(zzeta-1) );
C_star = C_star_bar;
I_star = I_star_bar;
R_star = R_star_bar;
PB_star = 1/(phi*R_star); 
Z = MC * (1-alp) * Yw / K;
RK =  ( Z + (1-delta) );
CH_star = C_star * ( 1/Terms )^(-muC_star);
IH_star = I_star  * ( 1/Terms )^(-muI_star);
UCR = (1-varrho)*(CR*(1-chi))^((1-varrho)*(1-sigma_c)-1)*(1-HR)^(varrho*(1-sigma_c));
YY=1;
HH=1;
WPWP=1;
CC=1;
II=1;
RR=1;
PiPi=1;
PiHPiH=1;
PiSPiS=1;
TT=1;
PFPF=1;
CHCH=1;
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ',1) = ' M_.param_names{iter,:} ';' ])
end




% The steady state valuues of the endogenous variables of the model are
% defined here.

NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for i = 1:NumberOfEndogenousVariables                         % Loop...
  %varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
  %eval(['ys(' int2str(i) ') = ' varname ';']);
  ys(i)=0;%    Get the steady state vZNue of this variable.
end          


