function y = int_new(xnodi, fnodi, x) 
%    Funzione che determina in un insieme di punti       
%    il valore del polinomio interpolante ottenuto        
%    dalla formula di Newton                               
%    Dati di input:                                        
%  
%    xnodi    vettore con i nodi dell'interpolazione         
%    fnodi    vettore con i valori nei nodi                  
%    x    vettore con i punti in cui si vuole            
%         calcolare il valore del polinomio              
%         interpolante                                  
%    Dati di output:                                       
%    y    vettore contenente i valori assunti            
%         dal polinomio interpolante                   
%                                                          
xnodi = xnodi(:); 
fnodi = fnodi(:); x = x(:);
n = length(xnodi);                                             
a = fnodi;                                                     
for i = 2:n                                                
   a(i:n) = (a(i:n)-a(i-1:n-1)) ...
      ./ (xnodi(i:n)-xnodi(1:n-i+1));  
end                                                        
y = a(n) * ones(size(x));                                       
for i = n-1:-1:1                                          
   y = y .* (x-xnodi(i)) + a(i);                            
end                                                        
