function [ys,check] = Main_steadystate(ys,exe)
global M_

%%
%% Param Values
%%
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(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;
%%
%% Exchange rates
%%
rer_H      =   1;
rer_U      =   1/rer_H;
%%
% Solver
 guess=abs(randn)*ones(4,1);
% guess=ones(6,1);
[res1 fval]     =   fsolve(@sfile, guess);
% while isreal(res1)==0
% guess=abs(randn)*ones(10,1);
% [res1]     =   fsolve(@sfile, guess);
% end
options = optimset('diagnostics','on', 'Display','off','FunValCheck','off','MaxFunEvals',2000, 'MaxIter', 3000,'DiffMaxChange',1,'TolFun',1e-10,'TolX',1e-10);

[res2 fval]    =   fsolve(@sfile, res1,options)
% fval
% res2

Y_H =res2(1);
Y_U =res2(2);
P_sc_H=res2(3);
P_sc_U=res2(4);
% Bs_H = res2(6);
% Bs_U = res2(7);

%%
%% Home
%%
P_sd_H = 1;

P_sm_H      =	rer_H*P_sc_H/P_sc_U;


P_si_H =  P_sc_H;%(   (1-dlt_i_H)*P_sd_H^(r_c / (r_c-1))  +   dlt_i_H*P_sm_H^(r_c /(r_c-1))  )^((r_c-1)/r_c);

r_k_H = P_si_H*(1/betta_H - (1-dlt));

w_H       = (   (0.5^(r_v/(r_v-1))    -   (1-dlt_l_H)*r_k_H^(r_v/(r_v-1)))/dlt_l_H  )^((r_v-1)/r_v);
mc_H =  (   (1-dlt_l_H)*r_k_H^(r_v/(r_v-1))    +   dlt_l_H*w_H^(r_v/(r_v-1))  )^((r_v-1)/r_v);

% mc_H      =  ((1-dlt_oy_H)*lmbd_sv_H^(r_o/(r_o-1))   +   dlt_oy_H*P_so_H^(r_o/(r_o-1))  )^((r_o-1)/r_o);

K_H  =   (1-dlt_l_H)*(mc_H/r_k_H)^(1/(1-r_v))*Y_H;
L_H  =  dlt_l_H*(mc_H/ w_H)^(1/(1-r_v))*Y_H;

theta_w_H=1;
psi_w_H	=	(1+theta_w_H)/theta_w_H;
w_f_H = w_H*(psi_w_H-1)/psi_w_H;   
C_H = (chi0*(P_sc_H /w_f_H)*(1-L_H)^(-chi))^(-1 /sig)*(1 /(1-kappa));
C_d_H =   (1-dlt_mc_H)*(P_sc_H)^(1/(1-r_c))*C_H;

I_H=dlt*K_H; 
I_d_H =  (1-dlt_i_H)*(P_si_H)^(1/(1-r_c))*I_H;

M_c_H   = dlt_mc_H*(P_sc_H/P_sm_H)^(1/(1-r_c))*C_H;
M_i_H   = dlt_i_H*(P_si_H/P_sm_H)^(1/(1-r_c))*I_H;
M_H = P_sm_H*(M_i_H + M_c_H);

G_d_H = g_H*Y_H;


%%
%%
P_sd_U = 1;
P_sm_U		=	rer_U*P_sc_U/P_sc_H;

P_si_U = P_sc_U;% (   (1-dlt_i_U)*P_sd_U^(r_c / (r_c-1))  +   dlt_i_U*P_sm_U^(r_c /(r_c-1))  )^((r_c-1)/r_c);

r_k_U = P_si_U*(1/betta_U - (1-dlt));

w_U       = (   (0.5^(r_v/(r_v-1))    -   (1-dlt_l_U)*r_k_U^(r_v/(r_v-1)))/dlt_l_U  )^((r_v-1)/r_v);
mc_U =  (   (1-dlt_l_U)*r_k_U^(r_v/(r_v-1))    +   dlt_l_U*w_U^(r_v/(r_v-1))  )^((r_v-1)/r_v);
% mc_U      =  ((1-dlt_oy_U)*lmbd_sv_U^(r_o/(r_o-1))   +   dlt_oy_U*P_so_U^(r_o/(r_o-1))  )^((r_o-1)/r_o);

K_U  =   (1-dlt_l_U)*(mc_U/r_k_U)^(1/(1-r_v))*Y_U;
L_U  =  dlt_l_U*(mc_U/ w_U)^(1/(1-r_v))*Y_U;


theta_w_U=1;
psi_w_U	=	(1+theta_w_U)/theta_w_U;
w_f_U = w_U*(psi_w_U-1)/psi_w_U;
C_U = (chi0*(P_sc_U /w_f_U)*(1-L_U)^(-chi))^(-1 /sig)*(1 /(1-kappa));
C_d_U =   (1-dlt_mc_U)*(P_sc_U)^(1/(1-r_c))*C_U;

I_U=dlt*K_U; 
I_d_U =  (1-dlt_i_U)*(P_si_U)^(1/(1-r_c))*I_U;

M_c_U   = dlt_mc_U*(P_sc_U/P_sm_U)^(1/(1-r_c))*C_U;
M_i_U   = dlt_i_U*(P_si_U/P_sm_U)^(1/(1-r_c))*I_U;
M_U = P_sm_U*(M_i_U + M_c_U);

G_d_U = g_U*Y_U;

%%
%%
X_H = M_U;
X_U = M_H;
NTs_H	=	(X_H	-	M_H);%	+	P_so_H*O_x;
NTs_U	=	(X_U	-	M_U );%	-	P_so_U*O_U;

r_s_H = 1/betta_H;   r_s_U = 1/betta_U;
P_b_H = betta_H;    P_b_U = betta_U;
Bs_H = NTs_H/(P_b_U - 1);
Bs_U = NTs_U/(P_b_H - 1);
% Bs_H = NTs_H/( 1 - r_s_U);
% Bs_U = NTs_U/( 1 - r_s_H);
%%
%%
%%
L_d_H = L_H;    L_d_U = L_U;    
omega_l_H = 1;
omega_H = 1;  
S1_w_H = L_H*psi_w_H*w_f_H/(w_H*(1-rho_w*betta_H*omega_H));
S2_w_H = L_H*(psi_w_H-1)*(1+tau_w)/(1 - rho_w*betta_H*omega_l_H);

omega_l_U = 1;
omega_U = 1;
S1_w_U = L_U*psi_w_U*w_f_U/(w_U*(1-rho_w*betta_U*omega_U));
S2_w_U = L_U*(psi_w_U-1)*(1+tau_w)/(1 - rho_w*betta_U*omega_l_U);

Lmbd_H = betta_H;    Lmbd_U = betta_U;  

pi_head_H =1; pi_head_U = 1;
Z_H=1;              Z_U=1;              
Z_c_H=1;            Z_c_U=1;            
Z_m_H=1;            Z_m_U=1;            
Z_i_H=1;            Z_i_U=1;            
Z_g_H=1;            Z_g_U=1;           

q_c_H = ((1-kappa)*C_H)^(-sig)/P_sc_H; 
q_c_U = ((1-kappa)*C_U)^(-sig)/P_sc_U; 
Q_k_H = betta_H*q_c_H*r_k_H/(1-betta_H*(1-dlt)); 
Q_k_U = betta_U*q_c_U*r_k_U/(1-betta_U*(1-dlt)); 

pi_d_H=1;
pi_l_H	=	pi_d_H^gama_p;
theta_p_H=1;  
psi_p_H	=	(1+theta_p_H)/theta_p_H;
S1_p_H = Y_H*psi_p_H*mc_H/(1-betta_H*rho_p*pi_d_H);
S2_p_H = Y_H*(psi_p_H - 1)*(1+tau_p)/(1-betta_H*rho_p*pi_l_H);
Y_d_H=Y_H;

pi_d_U=1;
pi_l_U	=	pi_d_U^gama_p;
theta_p_U=1;  
psi_p_U	=	(1+theta_p_U)/theta_p_U;

S1_p_U = Y_U*psi_p_U*mc_U/(1-betta_U*rho_p*pi_d_U);
S2_p_U = Y_U*(psi_p_U - 1)*(1+tau_p)/(1-betta_U*rho_p*pi_l_U);
Y_d_U=Y_U;


% Y_o = 1;
%%
%%
%%
% 
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 ';']);                %    Get the steady state value of this variable.
end