function I = simpson_rule(f, a, b, n)


    % Input validation
    if mod(n, 2) ~= 0
        error('n must be an even integer for Simpson''s Rule.');
    end

    % Calculate step size
    dx = (b - a) / n;

    % Initialize the sum with the first and last terms
    sum_f = feval(f, a) + feval(f, b);

    % Sum the odd indices (4 * f(x_i))
    for i = 1:2:n-1
        x_i = a + i * dx;
        sum_f = sum_f + 4 * feval(f, x_i); 
    end

    % Sum the even indices (2 * f(x_i))
    for i = 2:2:n-2
        x_i = a + i * dx; 
        sum_f = sum_f + 2 * feval(f, x_i); 
    end

    % Apply Simpson's rule formula
    I = (dx / 3) * sum_f;
end