function [ys, params, check] = PFE3_steadystate(ys, exo, M_, options_)
% Computes the steady state for the PFE model using a numerical solver
% Initialize the steady state computation check indicator
check = 0;
options = optimset('Display','iter','MaxIter',1000, 'Algorithm','levenberg-marquardt','TolFun',1e-8,'TolX',1e-8);  
 
% Extract parameter values
NumberOfParameters = M_.param_nbr;
params = M_.params; % Copying all parameter values into an array
 
% Assign parameter values to variables for easier access
for ii = 1:NumberOfParameters
    paramname = M_.param_names{ii};
    eval([paramname ' = M_.params(' num2str(ii) ');']);
end
 
 
 
% Define the steady state calculations
zeta=1;
a=1;
rhho = (sigma - 1) /sigma;
omega = 1 - n_r / n_w * (1 - gammad) ;
pi = pibar; 
pistar = ((1- theta * pi^(epsilon - 1)) / (1- theta))^(1 / (1-epsilon)); 
vp = ((1-theta) * pistar^(-epsilon)) / (1 - theta * pi^epsilon); 
r = 1 / betta - 1; 
i = pi * (1 + r) - 1; 
mc = ((epsilon - 1) / epsilon) * (1 - betta * theta * pi^epsilon) / (1 - betta * theta * pi^(epsilon-1)) * pistar; 
S = 0; 
S_prime = 0; 
q = 1 - S + S_prime; 
inv_k = delta; 
k_y = alppha / (q * (r + delta)) * mc * vp; 

%if k_y <= 0
 %   check = 1; % set failure indicator
  %  return;    % return without updating steady states
%end
 
L_y = (vp/(k_y^alppha))^(1 / (1 - alppha));  %0.208268818312174
 
%if l_y <= 0
   % check = 1; % set failure indicator
    %return;    % return without updating steady states
%end
 
inv_y = inv_k * k_y; 
c_y = 1 - inv_y; 
w = ((1 - alppha) * mc * vp) / L_y;  
mu = 1; 
r_rc = (1 + i) / (1 + i - gammad);
r_wc = ((1 - omega) * r_rc) / (1 + i - omega);
fr = frbar ;
tau = fr * nu * r_wc;
r_r = (1 + i) / (1 + i - mu * gammad);



 
 
    omegaw0 = 1;
    [omegaw, ~, exitflag] = fsolve(@(omegaw) omegaw - (omega + (1 - omega) * ((1 - tau + (nu / 8) * (mu * (omegaw - omega) * r_r) / ((1 + i) * omegaw - mu * omega)) * (1 / xi))^(1-v) * ((1 - gammad * betta^sigma * (1 + r)^(sigma - 1)) / (1 - betta^sigma * ((1 + r) * omegaw)^(sigma - 1)))^(1/(1 - sigma))), omegaw0, options);
    if exitflag <= 0
        check = 1; 
        return     
    end
 
 

 
 
r_w = (mu * (omegaw - omega) * r_r) / ((1 + i) * omegaw - mu * omega);
tau_w = tau - (r_w * nu / 8);
delta_r = 1 / (1 - gammad * (betta^sigma) * (1 + r)^(sigma - 1));
delta_w = 1 / (1 - (betta^sigma) * ((1 + r) * omegaw)^(sigma - 1));
lambda = v * delta_w^((rhho + 1) / rhho) * ((1 - v) / v * (1 / ((1 - tau_w) * w)))^(1 - v);

cw_y = ((1 - omega) * (n_w / n_r) * ((1 - tau_w) / xi)^(1 - v) * (delta_w / delta_r)^(sigma / (1 - sigma))) / ((betta * (1 + r) * omegaw)^sigma + (1 - omega) * (n_w / n_r) * ((1 - tau_w) / xi)^(1 - v) * (delta_w / delta_r)^(sigma / (1 - sigma)) - omega) * c_y;
cr_y = c_y - cw_y ;
 
M = ((betta * (1 + r) * omegaw)^sigma - omega) / ((1 - omega) * (n_w / n_r) * ((1 - tau_w) / xi)^(1 - v) * (delta_w / delta_r)^(sigma / (1 - sigma)));
 
KL = (alppha/(1-alppha)) * (w/(q*(r + delta)));

L = (n_w + xi * n_r) / ((M + (1 / (1 - tau_w))) * ((1-v) / (v * (1+M) * w)) * (KL^alppha - KL*delta)) ;


 
L_w = (L - xi * n_r + M * (1 - tau_w) * n_w) / (M * (1 - tau_w) +1) ;
 
L_r = n_r + (M / xi) * (1 - tau_w) * ((L - xi * n_r - n_w) / (M * (1 - tau_w) +1)) ;
 
 
 
 
d_r = n_r * xi * w;
h_r = ((1 + r) * d_r)/ (1 + r - gammad);

 

 
b_r = (1-omega) / (pi - gammad) * pi / (pi - omega) * (nu * w * L_w);
b_w = omega / (pi - omega) * nu * w * L_w;
a_c = (1+r)/r * (b_r - tau * w * L_w);
L_c = a_c /fr;
k= KL* L;
inv = delta * k;
y = L / L_y ;
f = (1 - mc * vp) * y;
c = c_y* y;
c_w = cw_y * y;
c_r = c - c_w;
x_r = c_r/v;
x_w= c_w/v;
d_w = n_w * (1 - tau_w) * w + f;
h_w = (d_w + (1 / (1 + r)) * (1 - omega / omegaw) * (n_w / n_r) * h_r) / (1 - omega / ((1 + r) * omegaw));
a_w = omega * ((1 - tau) * w * L_w + f - c_w) / ( 1 - omega * (1 + r));
a_r  = (c_r - mu * b_r - xi * w * L_r - (1 - omega) / omega * a_w) / r;
g2 = (lambda * y) / (1 - betta * theta * pi^(epsilon-1));
g1 = pistar * (epsilon - 1) / epsilon * g2;
 


 fprintf('pistar = %.9f\n', pistar);
fprintf('vp = %.9f\n', vp);
fprintf('r = %.9f\n', r);
fprintf('i = %.9f\n', i);
fprintf('mc = %.9f\n', mc);
fprintf('k_y = %.9f\n', k_y);
fprintf('L_y = %.9f\n', L_y);
fprintf('inv_y = %.9f\n', inv_y);
fprintf('c_y = %.9f\n', c_y);
fprintf('w = %.9f\n', w);
disp(['r_rc = ', num2str(r_rc)]);
disp(['r_wc = ', num2str(r_wc)]);
disp(['tau = ', num2str(tau)]);
disp(['r_r = ', num2str(r_r)]);
disp(['omegaw = ', num2str(omegaw)]);
disp(['tau_w = ', num2str(tau_w)]);
disp(['r_w = ', num2str(r_w)]);
disp(['delta_w = ', num2str(delta_w)]);
disp(['delta_r = ', num2str(delta_r)]);
disp(['lambda = ', num2str(lambda)]);
disp(['cw_y = ', num2str(cw_y)]);
disp(['cr_y = ', num2str(cr_y)]);
disp(['d_r = ', num2str(d_r)]);
disp(['h_r = ', num2str(h_r)]);
disp(['L_w = ', num2str(L_w)]);
disp(['c_w = ', num2str(c_w)]);
disp(['L_r = ', num2str(L_r)]);
disp(['L = ', num2str(L)]);
disp(['c_r = ', num2str(c_r)]);
disp(['h_w = ', num2str(h_w)]);
disp(['a_r = ', num2str(a_r)]);
disp(['a_w = ', num2str(a_w)]);
disp(['g1 = ', num2str(g1)]);
disp(['g2 = ', num2str(g2)]);
disp(['b_w = ', num2str(b_w)]);
disp(['b_r = ', num2str(b_r)]);
disp(['a_w = ', num2str(a_w)]);
disp(['c = ', num2str(c)]);
disp(['y = ', num2str(y)]);
disp(['f = ', num2str(f)]);
disp(['L_c = ', num2str(L_c)]);
disp(['d_r = ', num2str(d_r)]);
disp(['d_w = ', num2str(d_w)]);

 
 
%% Step 4: Update parameters and variables
params = NaN(M_.param_nbr, 1);
for iter = 1:M_.param_nbr
    eval(['params(' num2str(iter) ') = ' M_.param_names{iter} ';']);
end
 
for i = 1:M_.orig_endo_nbr
    eval(['ys(' int2str(i) ') = ' M_.endo_names{i} ';']);
end


