function [ x , F ] = F_integrale_Simpson( f , x0 , x_max , n )
  n = 2 * floor( n / 2 ) ;
  x = linspace( x0 , x_max , n+1 ) ;
  dx = x(2) - x(1) ;
  f_x = f( x ) ;
  p  = [ 5 8 -1 ] / 12 ;
  f1 = f_x(1:2:end-2) ;
  f2 = f_x(2:2:end-1) ;
  f3 = f_x(3:2:end  ) ;
  I_sx = p(1)*f1 + p(2)*f2 + p(3)*f3 ;
  I_dx = p(3)*f1 + p(2)*f2 + p(1)*f3 ;
  I = zeros( 1 , n ) ;
  I(1:2:end-1) = I_sx ;
  I(2:2:end  ) = I_dx ;
  F = [ 0 dx * cumsum( I ) ] ;
end