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

% Dominio
estremi_dominio = [ -6 6 ] ;
xv = 1 ;
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 e asintoti
f_dom = f( dominio ) ;
l = limite( f             , -10 ,  -1 , 1e-10 , false ) ; % x → -Inf
s = limite( f             ,  xv , -.1 , 1e-10 , true  ) ; % x →  xv-
d = limite( f             ,  xv , +.1 , 1e-10 , true  ) ; % x →  xv+
m = limite( @(x) f(x)/x   ,  10 ,  +1 , 1e-10 , false ) ; % x → +Inf
q = limite( @(x) f(x)-2*x ,  10 ,  +1 , 1e-10 , false ) ; % x → +Inf

% Plot
plot( dominio , f_dom , 'r' ) ;
hold on ;

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

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

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

xlabel('x') ;
ylabel('y') ;
ylim( [ 0 15 ] ) ;