TRES = trap    ( @fexample, 0, 10, 40 ) ;
SRES = simpson ( @fexample, 0, 10, 40 ) ;
ESA  = quad    ( @fexample, 0, 10 ) ;

fprintf( 1, 'trap %f\nsimp %f\nquad %f\nerr1 %f\nerr2 %f\n', ...
             TRES, SRES, ESA, TRES-ESA, SRES-ESA) ;

trapPlot( @fexample, 0, 10, 40) ;

H    = zeros(1,1) ;
TERR = zeros(1,1) ;
SERR = zeros(1,1) ;
k = 1 ;
for N=[10:10:60000]
  N
  H(k)    = (10-0)/N ;
  TERR(k) = trap   ( @fexample, 0, 10, N ) - ESA ;
  SERR(k) = simpson( @fexample, 0, 10, N ) - ESA ;
  k       = k+1 ;
end

axis square
loglog( H, abs(TERR), '-r', H, abs(SERR), '-b' ) ;
set(gca,'DataAspectRatio',[1 1 1],...
        'PlotBoxAspectRatio',[1 1 1])
    
