function Z_sol = solve_Z(tau_l, mc, mu, a, d_0, d_1, d_2, sigma, zeta, eta, delta_z, chi_e, phi_e, kappa, g_a_y_bar, chi_1, chi_2, g_y_bar)
   
    % Calcul des termes intermédiaires
    c_y = 1 - chi_1 * mu^chi_2 - g_y_bar;
    Z_y = ((1 + chi_e) * phi_e * (1 - mu) - kappa * g_a_y_bar) / (1 - delta_z);
    
    
    % Définition de l'équation
    equation_Z = @(Z) (1 - tau_l) * mc * (a * (1 - d_0 - d_1 * Z - d_2 * Z^2))^(sigma + 1) *  Z_y^(sigma + zeta)- eta * Z^(sigma + zeta) * c_y^zeta ;
    
    % Options de fsolve
    options = optimoptions('fsolve', 'Display', 'iter', 'TolFun', 1e-12, 'TolX', 1e-12, 'MaxIterations', 5000);
    
    % Tester plusieurs valeurs initiales
    Z0 = 500;
    
    % Résolution de l'équation pour Z
    [Z_sol, fval, exitflag] = fsolve(equation_Z, Z0, options);
    
    % Affichage de la solution pour Z
    fprintf('Solution pour Z: %.4f\n', Z_sol);
    
    % Vérification de la convergence
    if exitflag <= 0
        disp('fsolve n''a pas convergé vers une solution.');
    end
end