function mnhf_trapz_simpson_example(N)
%MNHF_TRAPZ_SIMPSON_EXAMPLE estimates the integral of cos(x) from x=0 to 
% x=10 using the built-in function trapz and a Simpson's rule solver. The 
% exact solution is sin(10).
%
% N - Number of grid points
%
% Usage  mnhf_trapz_simpson_example(11);
%        mnhf_trapz_simpson_example(101);
%        mnhf_trapz_simpson_example(1001);  
%        mnhf_trapz_simpson_example(10001);

x = linspace(0.0,10.0,N);
y = cos(x);
h = x(2)-x(1);
fprintf('h=%4.3f, Error (trapz)=%12.11f,',h,abs(sin(10)-h*trapz(y)))
fprintf(' Error (Simpson)=%12.11f\n',abs(sin(10)-mnhf_simpson(y,h)))

function area=mnhf_simpson(g,h)
%MNHF_SIMPSON finds the Simpson's Rule integral of discrete data that is 
% uniformly spaced on a grid.
%
% g - a vector containing the discrete data
% h - a scalar indicating the grid spacing

area = 0.0;       % Initialize sum.
n = length(g)-1;

if mod(n,2)==0  % Use Simpson's 1/3 rule.
 for ii=2:2:n
  area = area+4.0*g(ii);
 end
 for ii=3:2:n-1
  area = area+2.0*g(ii);
 end
 area = h/3.0*(g(1)+area+g(n+1));
else            % Use Simpson's 1/3 rule...
 for ii=2:2:n-3
  area = area+4.0*g(ii);
 end
 for ii=3:2:n-4
  area = area+2.0*g(ii);
 end
 area = h/3.0*(g(1)+area+g(n-2));
 % ... followed by Simpson's 3/8 rule. 
 area = area+0.375*h*(g(n-2)+3.0*g(n-1)+3.0*g(n)+g(n+1));
end