function s=Newton(x,y,x0,nn)
  % Newton插值，x與y為已知的插值點及其函數值
  % x0為要求的插值點的x坐標值。nn為Newton插值多項式的階數
  nx=length(x);
  ny=length(y);
  if nx~=ny
      warning('向量x與y的長度應該相同')
      return
  end
  m=length(x0);
  %按照公式，對要求的插值點x0的每個元素進行計算
  for i=1:m
      t=0.0;
      j=1;
      yy=y;
      kk=j;
      %求各級均差
      while(kk<=nn)
          kk=kk+1;
          for k=kk:nx
              yy(k)=(yy(k)-yy(kk-1))/(x(k)-x(kk-1));
          end
      end
      % 求差值結果
      t=yy(1);
      for k=2:nn
          u=1.0;
          jj=1;
          while(jj<k)
              u=u*(x0(i)-x(jj));
              jj=jj+1;
          end
          t=t+yy(k)*u;
      end
      s(i)=t;
  end