function [] = riemann(n)
% riemann: riemann sums demo
% riemann(n) calculates riemann sums using n steps

f = @(x) exp(-x^2); % monotone decreasing integrand

a = 0.0; % lower limit
b = 1.0; % upper limit

h = (b - a) / n; % rectangle widths

% Calculate upper and lower sums using fact
% that the integrand function is monotone decreasing

s = 0.0;
for i = n : -1 : 1
   s = s + f(a + i*h);
end

sumLower = s*h;
sumUpper = sumLower + h*(f(a) - f(b));

% Compare average of sums with exact value

average = 0.5*(sumLower + sumUpper);
exact = sqrt(pi)*erf(1.0)/2.0;

fprintf('Lower sum   = %.10f\n', sumLower);
fprintf('Upper sum   = %.10f\n', sumUpper);
fprintf('Average     = %.10f\n', average);
fprintf('Exact value = %.10f\n', exact);

end