function [I] = compmidpnt(f,a,b,N)
h = (b-a) ./ N;
x = linspace(a,b,N+1);
x_mid = (x(2:end)+x(1:end-1))/2;
I = sum(f(x_mid))*h;
end