function I = midpoint(a, b, f, N);

dx = (b - a) / N;
x = a + dx * (1 : N) - dx / 2;  % one point in the middle of each subinterval
y = f(x);
I = dx * sum(y);
