function I = trapIntegral(f, a, b, N)
% Compute the N point trapezoidal rule approximation of the function f
% on the interval [a, b]

x = linspace(a, b, N);
I = (b-a) * sum([f(x(1))/2, f(x(2:end-1)), f(x(end))/2])/(N-1);
end