X=poly(0,'x');

// points uniformément répartis

n=10;

// uniforme sur [-1,1]
// x = (-1:2/n:1);

// points de Chebyshev sur [-1,1]
k= (0:n)
x= cos( (2.*k+1)*%pi./(2*n+2) );

l = x;

// polynomes de base l_i

for i=1:n,
  l(i) = 1;
  for j=1:n,
    if (i<>j) then
      l(i)=l(i)*(X-x(j))/(x(i)-x(j));
    end   
  end
end

// calcul du polynome d'interpolation

p=0;

for i=1:n,
  p = p+sin(x(i))*l(i);
end


// affichage du polynome

p

// graphe du polynome

t = (-5:0.1:5);

plot2d(t,horner(p,t),1);
plot2d(t,sin(t),2);



