function [ Ih ] = trapezoidal_composite_quadrature( fun, a, b, M )
% TRAPEZOIDAL_COMPOSITE_QUADRATURE approximate the integral of a function in
% the interval [a,b] by means of the composite trapezoidal quadrature formula
%  [ Ih ] = trapezoidal_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

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}

Ih = H / 2 * sum( f_x_km1 + f_x_k );

return