function p=lagrange(x,y,punto);
% calcolo del polinomio di interpolazione di Lagrange
% relativo ai vettori di punti (x,y), in punto

%cerco se punto e' uno dei valori presenti nel vettore x
k=find(x==punto);
%se lo e' non devo fare altro che assegnare a p il valore y corrispondente
if length(k)~=0
   p=y(k(1));
%se non lo e' utilizzo Lagrange
else
   n=length(x)-1;
   p=0;
   z=prod(punto-x); %z=produttoria(x-x(k))
   for i=0:n
      r=z;
      for j=0:n
         if i~=j
            r=r/(x(i+1)-x(j+1)); %r=(x-x(j))/(x(k)-x(j))
         endif;
      endfor;
      p=p+y(i+1)*r/(punto-x(i+1)); %eseguo la sommatoria per il caso generale di Lagrange
   endfor;
endif;
%fine programma
