function integral = int_trapezoid ( f, a, b, N)
% integrates f(x) in the a to b range with trapezoid method
% over N points

if( N ~= floor(N) ) error('N must be integer'); end
if( N < 2 ) error('N>=2 is required'); end

x=linspace(a,b,N);
h=(b-a)/(N-1);

integral = (f(a)+f(b))/2; % end points are special

for i=2:N-1
	xi=a+h*(i-1);
	integral = integral+f(xi);
end

integral=integral*h;

end
