function [M,a]=NDD(x,f)
%function [M,a]=NDD(x,f)
%
%  Inputs:
%    x - vector of x values
%    f - vector of funtion values corresponding to the x values
%  Outputs:
%    M - matrix of divided differences
%    a - coeficients for divided difference polynomial
%
%  See also: NDDPolyEval.m
M=[x, f];
n=length(x);
a=f(1);
v=zeros(size(x));
f2=f;
for j=1:(n-1)
    for k=1:j
        v(k)=0;
    end
    for k=(j+1):n
        v(k)=(f2(k)-f2(k-1))/(x(k)-x(k-j));
    end
    M=[M, v];
    a=[a; v(j+1)];
    f2=v;
end