function area = trapezoid( f, a, b, n )
%TRAPEZOID Integrate f from a to b via the trapezoid rule on n intervals.


% Set up x and f(x) values. Don't assume that f vectorizes.

x = linspace( a, b, n+1 );

y = [];
for i = 1 : n+1
    y(i) = f( x(i) );
end


% Calculate the area under the curve

deltaX = ( b - a ) / n;
area = sum( deltaX * ( y(1:n) + y(2:n+1) ) / 2 );


% Plot the trapezoids by listing (x,y) coordinates for each corner of
% each trapezoid and then plotting a line that connects all those
% coordinates.

trapX = [];
trapY = [];

for i = 1 : n
    trapX = [ trapX x(i) x(i) x(i+1) x(i+1) ];
    trapY = [ trapY 0 y(i) y(i+1) 0 ];
end

clf;
plot( trapX, trapY, 'r' );
hold on;


% Plot the function using a more dense set of samples than the trapezoids
% probably used. Once again, don't assume the function vectorizes.

nSamples = 200;

funcX = linspace( a, b, nSamples );
funcY = [];

for i = 1 : nSamples
    funcY(i) = f( funcX(i) );
end

plot( funcX, funcY, 'b' );

end

