function I = midpt(f, a, b, n)
if nargin < 4, n = 100; end
h = (b - a) / n;
x = a:h:b;
s = 0;
for i = 1:(n-1)
    s = s + f(x(i) + h/2);
end
I = h * s;
end
