n=10; % il grado
theta=linspace(0,pi,n+1); % gli angoli
x=cos(theta); % i punti di Chebyshev-Lobatto
f=@(x)sin(2*pi*x) % la funzione da interpolare 
y=f(x); % i valori da interpolare
eps=10^(-10);
t=linspace(-1+eps,1-eps,1017); % modello discreto di [-1,1]
%
% prima il polinomio d'interpolazione usando polyfit
%
tic
p=polyfit(x,y,n);
vp=polyval(p,t);
toc
plot(t,f(t),t,vp)
%
% adesso la forma baricentrica
%
tic
s1=y(1)./(t-x(1))/2;
s2=1./(t-x(1))/2;
for i=1:(n-1)
    s1=s1+(-1)^i*y(i+1)./(t-x(i+1));
    s2=s2+(-1)^i./(t-x(i+1));
end
s1=s1+(-1)^n*y(n+1)./(t-x(n+1))/2;
s2=s2+(-1)^n*1./(t-x(n+1))/2; 
vb=s1./s2;
toc
figure
plot(t,f(t),t,vb)