function T = traprule(f,a,b,n)
% TRAPRULE uses n intervals on [a,b] and does the trapezoid rule
% try:
%   >> traprule(@(x) sin(x), 0, 2, 4)
% and compare (use "format long" for more digits)
%   >> simprule(@(x) sin(x), 0, 2, 4)
%   >> quad(@(x) sin(x), 0, 2)          % Matlab built-in
%   >> - cos(2) + cos(0)                % exact

dx = (b - a) / n;    % user provides  a, b, n
T = f(a);
for i = 1:n-1
  x = a + i * dx;
  T = T + 2 * f(x);
end
T = T + f(b);
T = ((b - a) / (2 * n)) * T;
