function c = newtonkoeff(x,y)



%berechnet zu den daten (x(1),y(1)), ..., (x(n),y(n)) die
%koeffizienten des interpolationspolynoms von hoechstgrad n-1 in der form
%von Newton p(x)=((((((x-x(n-1))*c(n)+c(n-1))*(x-x(n-2)+...)*(x(1))+c(1)
n = length(x);
if length(y)~=n
  error ('newtonkoeff: length(x) not equal length(y)');
end

divdiff=zeros(n,n);
c=divdiff(1,:);
divdiff=zeros(n,n);
divdiff(:,1)=y;
for j=2:n
   divdiff(1:n-j+1,j)=(divdiff(2:n-j+2,j-1));
        %...divdiff(1:n-j+1,j-1))./(x(j:n)-x(1:n-j+1));
   end
c=divdiff(1,:);