
function R_solution1 = solve_for_R_tau_y(n, tau_e, tau_k, b_y)

g_y_bar   =0;
delta_k   = 0.1; 
alpha     = 1/3;
chi_1     = 0.1;
chi_2     = 2.6; 
phi_e     = 0.45;
beta      = 0.985; 

  
    % Calcul du k_y en fonction de R1
    k_y = @(R1) alpha * (1 - phi_e*(1-((tau_e*phi_e/(chi_1*chi_2))^(1/(chi_2-1))))*tau_e - chi_1*((tau_e*phi_e/(chi_1*chi_2))^(1/(chi_2-1)))^chi_2) *(1 - tau_k)/(R1-1+delta_k);
    
      % Calcul du c_y en fonction de R1
    c_y = @(R1)1-chi_1*((tau_e*phi_e/(chi_1*chi_2))^(1/(chi_2-1)))^chi_2 -g_y_bar - (n+delta_k)*k_y(R1);
    
    % Définir l'équation pour R1
    equation_R = @(R1) c_y(R1) * (R1-1/beta) - n * (1 - beta) / beta * R1*(b_y + k_y(R1));
    
    % Guess initial pour R
    R_initial_guess = 2;
    
    % Résoudre pour R en utilisant fsolve
    options = optimoptions('fsolve', 'Display', 'off', 'TolFun', 1e-9, 'TolX', 1e-9);
    [R_solution1, fval, exitflag] = fsolve(equation_R, R_initial_guess, options);
    
    % Vérification du succès de la résolution
    if exitflag <= 0
        error('FSOLVE did not converge to a solution.');
    end
    
    % Diagnostic de sortie
    fprintf('FSOLVE exit flag: %d, fval: %f\n', exitflag, fval);
end
