simpsons rule using tolerance in matlab

main.m
function [approximation, n] = simpsons_rule_tolerance(f, a, b, tol)
% SIMPSONS_RULE_TOLERANCE Approximate definite integral of a function using Simpson's rule with tolerance.
%   approximation = SIMPSONS_RULE_TOLERANCE(f, a, b, tol) approximates the definite integral of the
%   function 'f' over the interval [a,b] using Simpson's rule with a relative tolerance 'tol'.
%   The function returns the approximation and the number of iterations 'n' used to reach the tolerance.

    % Step 0: set initial values
    n = 1; % number of intervals
    h = (b-a)/2; % interval width
    x = [a, (a+b)/2, b]; % nodes
    y = f(x); % function values at nodes
    approx_old = h/3*(y(1)+4*y(2)+y(3)); % compute initial approximation

    while true
        % Step 1: divide intervals and compute new approximation
        n = 2*n; % double the number of intervals
        h = h/2; % halve the interval width
        x_new = a:h:b; % nodes for new intervals
        y_new = f(x_new); % function values at new nodes
        approx_new = h/3*(y_new(1) + 4*sum(y_new(2:2:end-1)) + 2*sum(y_new(3:2:end-2)) + y_new(end)); % compute new approximation

        % Step 2: check tolerance
        if abs(approx_new - approx_old) < tol*abs(approx_new)
            approximation = approx_new; % approximation within tolerance
            return;
        end
        
        % Step 3: update values and approximation
        x = x_new; % update nodes
        y = y_new; % update function values at nodes
        approx_old = approx_new; % update approximation
    end
end
1550 chars
34 lines

main.m
f = @(x) sin(x); % define function
a = 0; b = pi/2; % interval endpoints
tol = 1e-6; % relative tolerance
[approximation, n] = simpsons_rule_tolerance(f, a, b, tol); % use Simpsons rule with tolerance
fprintf('Approximation: %.8f\n', approximation); % print approximation
fprintf('Intervals: %d\n', n); % print number of intervals used
336 chars
7 lines

main.m
Approximation: 0.99999923
Intervals: 16
40 chars
3 lines