function S = compSimp(f,a,b,n)
% Composite trapezoid rule for approximating the definite integral of a
% function.
%
% INPUT: f (function to be integrated)
%        a (left endpoint)
%        b (right endpoint)
%        n (even number of subintervals)
%
% OUTPUT: T (approximate integral value)
%
if mod(n,2) ~= 0
    fprintf('n must be an even positive number');
end

h = (b-a)/n; % Computing stepsize
x = linspace(a,b,n+1); % Vector of nodes
y = f(x); 
S = y(1)+y(end);
S = S + 4*sum(y(2:2:end-1));
S = S + 2*sum(y(3:2:end-2));
S = h*S/3;