function [final_result, n] = Simpson_converge_6_decimal(f, a, b)
    % Analytical result (log(2)) for comparison
    analytical_result = log(2);
    
    % Initialize variables
    n = 1; % Start with 1 interval
    final_result = 0; % To store the final result
    accuracy_threshold = 1e-6; % Desired accuracy for 6 decimal places
    
    % Loop until the result matches the analytical result to 6 decimal places
    while true
        % Calculate the integral using Simpson's rule
        delta_x = (b - a) / n;
        f_a = f(a);
        f_b = f(b);
        sum_f_i_odd = 0; % For 4 * sum(f_i) where i is odd
        sum_f_i_even = 0; % For 2 * sum(f_i) where i is even

        % Calculate function values for all intervals
        for i = 1:n-1
            x_i = a + i * delta_x; % Midpoint for each sub-interval
            f_i = f(x_i);
            if mod(i, 2) == 1 % odd index
                sum_f_i_odd = sum_f_i_odd + f_i;
            else % even index
                sum_f_i_even = sum_f_i_even + f_i;
            end
        end

        % Combine the sums into the final approximation
        current_result = (delta_x / 3) * (f_a + 4 * sum_f_i_odd + 2 * sum_f_i_even + f_b);
        
        % Check if current result matches analytical result to 6 decimal places
        if round(current_result, 6) == round(analytical_result, 6)
            final_result = current_result;
            break; % stop when the result is accurate to 6 decimal places
        end

        % Update for next iteration
        n = n + 1; % increase the number of intervals
    end

    % Display the result and the number of intervals needed
    fprintf('Final result: %.6f\n', final_result);
    fprintf('Number of intervals needed: %d\n', n);
end