function yd = NewtonDiff(x,y,xd)
% Differentiation algorithm based on Newton interpolation
% x,y - row-vectors of (n+1) data values (x,y)
% xd - a row-vector of x-values, where numerical derivatives are to be found 
% yd - a row-vector of the numerical derivatives at the points of (xd)
% the numerical derivatives have order of O(h^n) in this algorithm

n = length(x) - 1; % the degree of interpolation polynomial
nd = length(xd); % the number of numerical derivatives to be found

D = ones(n,n); % the matrix for Newton divided differences 
                    % has n rows for divided differences of the corresponding order
                    % has n columns corresponding to data values x(1),x(2),...,x(n)
                   
for k = 1 : n 
    D(k,1) = (y(k+1)-y(k))/(x(k+1)-x(k)); % first divided difference
end
                   
for k = 2 : n % start the outer loop through the k-th order divided differences
    
    for kk = 1 : (n-k+1) % start the inner loop through the data values 
                        % the matrix of divided differences is triangular!
                        
        D(kk,k) = (D(kk+1,k-1) - D(kk,k-1))/(x(kk+k)-x(kk)); % see the Newton divided differences
    end % end of the inner loop
        
end % the end of the outer loop

% define the (n+1) coefficients of the Newton interpolating polynomial
a(1) = y(1); 

for k = 2 : (n+1)
    a(k) = D(1,k-1); % the first row of the matrix of divided differences corresponds to the point x(1)
end

% compute the values of numerical derivatives from the derivative of the 
% Newton interpolating polynomial at each x-value of xd

% Note: the algorithm uses double nested multiplication for polynomial evaluation
% Note: the algorithm also uses the matrix capacities of MATLAB

ydd = a(n+1)*ones(n,nd); 
% ydd - the matrix of n rows for n nested multiplications and nd columns for nd points of xd
% each element of yd is initialized as the coefficient of the highest degree of the polynomial

for j = 1 : (n-1) % the loops through rows of yd 

	for k = 1 : (n-j) % the loop to compute each nested multiplication for all values of xd in parallel
	    ydd(j,:) = a(n+1-k)+ydd(j,:).*(xd-x(n+1-k)); 
	end
end

yd = ydd(n,:); 
% initialization of the vector yd obtained in a final nested multiplication

for k = 1 : (n-1)
    yd = ydd(n-k,:) + yd.*(xd-x(n-k)); 
end





