f = @(x) sin( 7/2 * x ) + exp( x ) - 1;   a = 0;   b = 1; 
I_ex = exp( 1 ) - 2 + 2 / 7 * ( 1 - cos( 7 / 2 ) )
%   I_ex =
%        1.2716
M = 1;
Imp = midpoint_composite_quadrature( f, a, b, M )
%    Imp =
%       1.6327
It = trapezoidal_composite_quadrature( f, a, b, M )
%    It =
%        0.6837
Is = simpson_composite_quadrature( f, a, b, M )
%    Is =
%        1.3164

M = 10;
Imp_c = midpoint_composite_quadrature( f, a, b, M )
%    Imp_c =
%        1.2737
It_c = trapezoidal_composite_quadrature( f, a, b, M )
%    It_c =
%        1.2673
Is_c = simpson_composite_quadrature( f, a, b, M )
%    Is_c =
%        1.2716

M_vect = 2.^[2:1:7];
H_vect = ( b - a ) ./ M_vect; 
Err_mp_c = [ ];  Err_t_c = [ ];   Err_s_c = [ ];
for M = M_vect
    Imp_c = midpoint_composite_quadrature( f, a, b, M );
    Err_mp_c = [ Err_mp_c, abs( Imp_c - I_ex ) ];    
    It_c = trapezoidal_composite_quadrature( f, a, b, M );
    Err_t_c = [ Err_t_c, abs( It_c - I_ex ) ];    
    Is_c = simpson_composite_quadrature( f, a, b, M );
    Err_s_c = [ Err_s_c, abs( Is_c - I_ex ) ];
end
% Plot of the errors and the reference curves H^2 and H^4
loglog( H_vect, Err_mp_c, '-ob', H_vect, Err_t_c, '-xr', H_vect, Err_s_c, '-sm',...
            H_vect, H_vect.^2, '--k', H_vect, H_vect.^4, '-.k' );
legend( 'E_{mp}(f)', 'E_{t}(f)', 'E_{s}(f)', 'H^2', 'H^4');

conv_order_mp_c = log( Err_mp_c( end ) / Err_mp_c( end - 1 ) ) / ...
                         log( H_vect( end ) / H_vect( end - 1 ) )
%   conv_order_mp_c =
%       2.0001
conv_order_t_c = log( Err_t_c( end ) / Err_t_c( end - 1 ) ) / ...
                         log( H_vect( end ) / H_vect( end - 1 ) )
%   conv_order_t_c =
%       2.0001
conv_order_s_c = log( Err_s_c( end ) / Err_s_c( end - 1 ) ) / ...
                         log( H_vect( end ) / H_vect( end - 1 ) )
%   conv_order_s_c =
%       4.0001

f = @(x,d) x.^d; a=0;  b=1;
Table = [ ];
for d = [ 0 : 4 ]    
    I_ex = 1 / ( d + 1 );  % exact integral
    Imp = midpoint_composite_quadrature( @(x)f(x,d), a, b, 1 );    
    It = trapezoidal_composite_quadrature( @(x)f(x,d), a, b, 1 );
    Is = simpson_composite_quadrature( @(x)f(x,d), a, b, 1 );
    err_mp = I_ex - Imp;     err_t = I_ex - It;        err_s = I_ex - Is;
    Table = [ Table; [ d, err_mp, err_t, err_s ] ];
end
disp('         d, err_{mp},    err_t,    err_s'); disp( Table );

%            d, err_{mp},    err_t,    err_s
%            0         0         0         0
%       1.0000         0         0         0
%       2.0000    0.0833   -0.1667         0
%       3.0000    0.1250   -0.2500         0
%       4.0000    0.1375   -0.3000   -0.0083