function [c,w]=divdiff(x,y)

% Divided difference of data y with respect to the knots in x
%
% x,y: vectors of equal length
%
% Output:
% c: the divided difference 
% w: the vector of coefficients for the Newton polynomial
%
% repeated x-values are allowed if x is increasing
% then the data values in y are taken as function,
% derivative and higher derivative values, respectively.

[x,i]=sort(x);
y=y(i);
l=length(y);
if (l ~= length(x))
   error('x and y must have the same dimension.')
end

%use the spline toolbox to find multiplicities
m=myknt2mlt(x); % see description for knt2mlt
index=find(m>0);
ys=y; %save for derivative values
y(index)=y(index-m(index)); % fill function values

for k=1:l-1
   % for differences when consecutive x-values are different
   xdiff=(x(k+1:l)-x(1:l-k));
   ydiff=diff(y(k:l));
   index=find(m(k+1:l)<k);
   y(k+index)=ydiff(index)./xdiff(index);
   % for other values insert the proper derivative
   index=find(m>k-1);
   y(index)=ys(index-m(index)+k);
end

w=y;
c=y(end);


function knmult=myknt2mlt(t)
%
% Berechnung der Vielfachheit eines Knotens in einem 
% aufsteigend sortierten Knotenvektor
% knmult(i) ist die genaue Anzahl der Einträge t(k)=t(i)
% mit k<i, also 
%   knmult(i)=0 für das erste Auftreten von t(i) (also t(i-1)<t(i)),
%            =1 für das zweite Auftreten (also t(i-2) < t(i-1)=t(i),
%   etc.
lt=length(t);
difft = diff(t); % ergibt Null für gleiche Knoten

index = zeros([1,lt]); 
index(2:end) = (difft==0);
knmult = cumsum(index);

zz = find(diff(index)<0); % hier springt man zum nächsten Knoten
if isempty(zz) return, end % es lag nur ein Knoten der Vielfachheit lt vor

z = zeros([1,lt]);
pt = knmult(zz);
pt(2:end) = diff(pt);
z(zz+1) = pt;
knmult = knmult - cumsum(z);