// Interpolation de Lagrange
// Phenomene de Runge
//
n=evstr(x_dialog('nombre de points d''interpolation','10'));
//deff('y=f(x)','y=exp(-x)');  // un cas qui marche
deff('y=f(x)','y=ones(x)./(ones(x)+x.^2)')//un cas qui ne marche pas 
//
l=evstr(x_dialog('taille de l intervalle centre en zero ?','10')); 
a=-l/2;h=l/n;
//
// polynome d'interpolation calculé par l'algorithme des différences divisées
//
for j=1:n+1
 d(j)=f(a+(j-1)*h);
end;
b(1)=d(1);
for k=1:n
  for j=1:n+1-k
    d(j)=(d(j+1)-d(j))/(k*h);
  end
  b(k+1)=d(1);  
end
deff('y=p(x)',['y=b(n+1);for k=n:-1:1';
    'y=b(k)+(x-a-(k-1)*h)*y;end']);
//      
x=-l/2:l/100:l/2;
clf()
fplot2d(x,p,1);
fplot2d(x,f,2);
fplot2d(-l/2:h:l/2,f,-2);
legends(['polynome','fonction'],[1,2])
//
