%% hump_simpson applies composite Simpson rules to the hump() function.
%
a = 0.0;
b = 2.0;

fprintf ( 1, '\n' );
fprintf ( 1, '     n         q                  e\n');
fprintf ( 1, '\n' );

for nlog = 0 : 10
  n = 2 ^ nlog;
  x = linspace ( a, b, 2 * n + 1 );

  q =              hump(x(1) ) ...
     + 4.0 * sum ( hump(x(2:2:2*n)) ) ...
     + 2.0 * sum ( hump(x(3:2:2*n-1)) ) ...
     +             hump(x(2*n+1));

  q = ( b - a ) * q / 6.0 / n;
  e = hump_int ( a, b ) - q;
  fprintf ( 1, '  %4d  %20.16g  %10.4g\n', n, q, e );
end

