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

dx = (b - a) / N;
x = a + dx * (1 : N-1);
y = f(x);
I = dx / 2 * (f(a) + 2*sum(y) + f(b));
