function [ Ih ] = simpson_composite_quadrature( fun, a, b, M )
% SIMPSON_COMPOSITE_QUADRATURE approximate the integral of a function in
% the interval [a,b] by means of the composite Simpson quadrature formula
%  [ Ih ] = simpson_composite_quadrature( fun, a, b, M )
%  Inputs: fun = function handle, 
%          a,b = extrema of the interval [a,b]
%          M = number of subintervals of [a,b] of the same size, M>=1 
%              (the case M=1 corresponds to the simple formula)
%  Output: Ih = approximate value of the integral
%

H = ( b - a ) / M;

x_k = linspace( a, b, M + 1 ); % M+1 nodes
x_bar_k = ( x_k( 1 : end - 1 ) + x_k( 2 : end ) ) / 2; % M coordinates

f_x_km1 = fun( x_k( 1 : end - 1 ) ); % M values, in x_{k-1}
f_x_k = fun( x_k( 2 : end ) ); % M values, in x_{k}
f_x_bar_k = fun( x_bar_k ); % M values, in x_{k_bar}

Ih = H / 6 * sum( f_x_km1 + 4 * f_x_bar_k + f_x_k );

return