function pol=newton_pol(z,x,y)
  D=divdiff(x,y);
  n=size(D,2);
  pol=D(1,2);
  for k=3:n
  zz=1;
    for i=1:k-2
      zz=zz.*(z-x(i));
    end
    pol=pol+zz.*D(1,k);
  end
  