n=input('Grado del polinomio= ');
a=-5; b=5;
fun='1./(1+x.^2)';
%% Nodi equispaziati
x=linspace(a,b,n+1);
y=eval(fun);
p=diffdiv(x,y);
%% Nodi di Chebyshev
x=-cos([0:n]*pi/n);
x=(a+b)/2+(b-a)/2*x;
y=eval(fun);
chp=diffdiv(x,y);
%% Grafico della funzione e dei due polinomi interpolatori
x=linspace(-5,5);
fx=eval(fun);
px=polyval(p,x);
chpx=polyval(chp,x);
plot(x,fx,x,px,x,chpx)
legend('Funzione','Nodi equispaziati','Nodi Chebyshev')