% Definizione funzione
g = @(x) (x.^4+3*x.^3+2*x.^2-x-4) ./ (x.^3-1) + x.*tanh(x) ;

% Dominio di interesse
estremi_dominio = [ -6 6 ] ; % estremi
xv = 1 ;                     % valore di x che annulla il denominatore

% Definizione derivata approssimata per determinare gli intervalli di partenza
h_approx     = 1e-3 ;
dg_dx_approx = @(x) ( g(x+h_approx)-g(x-h_approx) ) / (2*h_approx) ;

% Determinazione degli intervalli di partenza con metodo del cambio del segno
N = 1000 ;

%  zeri a sx dell'asintoto verticale xv = 1
ab_sx = [ estremi_dominio(1) xv ] ; % intervallo
 z_sx = metodo_cambio_segno( dg_dx_approx , ab_sx , N ) ;

%  zeri a dx dell'asintoto verticale xv = 1
ab_dx = [ xv estremi_dominio(2) ] ; % intervallo
 z_dx = metodo_cambio_segno( dg_dx_approx , ab_dx , N ) ;

%  concatenazione zeri a sinistra e a destra dell'asintoto xv = 1
z = [ z_sx z_dx ] ;

% Derivata mediante limite
 d_dx = @(f,x0) limite( @(h) ( f(x0+h)-f(x0-h) ) / (2*h) , 0 , .1 , 1e-6 , true ) ;
dg_dx = @(x) d_dx(g,x) ;

% Zeri della derivata mediante metodo delle secanti applicato a dg_dx
n_zeri = size(z, 2) ;
zeri   = zeros(1, n_zeri) ;
for i = 1 : n_zeri
  zeri(i) = metodo_secanti(dg_dx, z(:,i), 1e-6) ;
end

% Plot funzione g
N  = 1000 ;
dominio_sx = linspace( estremi_dominio(1) , xv , N ) ; % [ -6 , 1 ]
dominio_dx = linspace( xv , estremi_dominio(2) , N ) ; % [  1 , 6 ]
dominio    = [ dominio_sx , dominio_dx ] ;

%  funzione
f_dom = g( dominio ) ;
plot( dominio , f_dom , 'r' ) ;
hold on ;

%  asintoto orizzontale x → -Inf
x = [ estremi_dominio(1) , -2 ] ;
y = 3 + 0*x ;
plot( x , y , '--k' ) ;

%  asintoto obliquo x → +Inf
x = [ 2 , estremi_dominio(2) ] ;
y = 2*x + 3 ;
plot( x , y , '--k' ) ;

%  asintoto verticale x → xv
y = [ -1000 1000 ] ;
plot( [ xv xv ] , y , '--k' ) ;

% Assi e limitazione del grafico lungo y
%xlabel('x') ;
%ylabel('y') ;
ylim( [ 0 15 ] ) ;

% Calcolo dei valori di massimo/minimo di g
g_minmax = g( zeri ) ;

% Plot con marker circolari blu:
plot( zeri , g_minmax , 'ob' ) ;

% Definizione derivata seconda
ddg_dx = @(x) d_dx(dg_dx,x) ;

% Plot dell'approssimante quadratico nei punti di estremo
x_quadratico = linspace(-1, 1, 400) ;
for i = 1 : n_zeri
  ddg_i = ddg_dx( zeri(i) ) ;
    x_i = x_quadratico + zeri(i) ;
    y_i = g_minmax(i) + (ddg_i/2) * x_quadratico .^ 2 ;
    plot( x_i , y_i , '-.b' ) ;
end