function [int, err] = trap_demo(a, b, m)
   h = (b - a)/m;
   x = a:h:b;
   y = sin(x);
   q = cumtrapz(x, y);
   int = q(length(q));
   exact = -cos(b) + cos(a); % exact = int_a^b sin(x) dx
   err = abs(exact - int);
end