% The function y allows you to calculate the value of the % interpolating polynomial, given some starting points in % arbitrary values of the variable
% independent.
 
%It accepts INPUT% in three vectors:
% x: whose components constitute the abscissas of the %points interpolating;
% y whose components are the ordinates of the %interpolating points;
% t: tab vector, whose components are the values on which
% we want to calculate the polynomial interpolating the %points (x (i), f (i))
% and a scalar:
% TOLL: tolerance within which must be the absolute value % of
% difference of the abscissa (or ordinate) of the points % of interpolation.
% In OUTPUT instead provides the vector y whose  %components %are the calculated values of the polynomial 
% in the components of the vector t tab.

function y= Newton(x, f, t, TOLL)
n= length(x); m=length(f); k =length(t); y=1:1:k;
if n~= m % length control vectors x and f
    error('error:interpolation points unspecified');
end
   A=zeros(n,n); % initialize matrix designed to contain %the divided differences
   for i=1:1:n
       A(i,i)=f(i);
   end
   for i=2:1:n
       for j=i-1 :-1:1
           if abs(x(i)-x(j))<= TOLL && abs(f(i)-f(j))<= TOLL
               error('coincident points');
           elseif abs(x(i)-x(j))<= TOLL
               error('abscissas coincident')
           end
           A(i,j)= (A(i,j+1) - A(i-1,j))/(x(i)-x(j)); 
           % the matrix A at the end of the cycle will be %inferiorly triangular
           % and is such that for i> j  %A(i,j)=f[x_j,x_(j+1),....,x_i]
           %with f [, ..] is the divided difference
       end
   end
   for h=1:1:k % calculation of the interpolating %polynomial in the components of t
       p=A(n,1);
       for i=n-1:-1:1
           p=(t(h) - x(i))*p + A(i,1);
       end
       y(h)=p;
   end
end