function [ Ih ] = I2c( fun, a, b, N )
%  I2c approximate the integral of a function in the interval 
%  [a,b] by means of the composite quadrature formula I_2^c 
%  defined in exercise 2, question 3. 
% 
%  [ Ih ] = I2c( fun, a, b, N )
%  Inputs: fun = function handle, 
%          a,b = extrema of the interval [a,b]
%          N = number of subintervals of [a,b] of the same size, N>=1 
%              (the case N=1 corresponds to the simple formula)
%  Output: Ih = approximate value of the integral
%

H = ( b - a ) / N;

x_k = linspace( a, b, N + 1 ); 
x_bar_k = ( x_k( 1 : end - 1 ) + x_k( 2 : end ) ) / 2;

f_x_bar_k = fun( x_bar_k );
f_x_m1 = fun( x_k( 1 : end - 1 ) );
f_x_p1 = fun( x_k( 2 : end ) );

Ih = H * sum( f_x_p1 + 4*f_x_bar_k + f_x_m1 ) / 6;

return