% Input: f, function handle alla funzione f
%        ab, vettore [ a b ] dell'intervallo entro al quale lavorare
% Output: xd, matrice 2 x N, N = numero di punti di discontinuità
%             Ogni colonna di xd è un intervallo che racchiude un punto di
%             discontinuità
function xd = trova_discontinuita( f , ab )
  dx  = 1e-3 ;
  x   = ab(1) : dx : ab(2) ;
  f_x = f( x ) ;
  ddf = ( f_x(1:end-2) - 2*f_x(2:end-1) + f_x(3:end) ) / (dx^2) ;
  M   = 1e3 ;
  ix  = find( abs(ddf) > M ) ;
  xd  = x(ix+1) ;
  xd  = reshape( xd , 2 , length(xd) / 2 ) ;
end