function S = Simpson(f, a, b, n)
    % integrate: Approximates the integral of a function f from a to b
    % using a specific numerical integration method.
    %
    % Inputs:
    %   f - function handle representing the function to integrate
    %   a - lower limit of integration
    %   b - upper limit of integration
    %   n - number of sub-intervals
    %
    % Output:
    %   I - approximated integral value

    % Calculate Δx
    delta_x = (b - a) / n;
    
    % Initialize the function values
    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
    S = (delta_x / 3) * (f_a + 4 * sum_f_i_odd + 2 * sum_f_i_even + f_b);
end