function [sum] = trap(f, a, b, n )

%trap trapezoidal integration
%   f is the function to integrate
%   a, b are the lower and upper integration limits
%   n is the number of steps

h = (b - a) / n;
sum = 0.5*(f(a) + f(b));
for k = 1 : n-1
   sum = sum + f(a + k*h);
end
sum = h*sum;
end