function[y0,D]=newton_poly(X,Y,x0)%Input-X插值节点的自变量取值%-Y插值节点的函数值%-x0是插值点%Output-y0是牛顿插值多项式在x0的值%-D是均差表n=length(X);m=length(Y);%ifn~=m%error('ThelengthsofXandYmustbeequal!');%return%endD=zeros(n,n);D(:,1)=Y';%Formthedivided-differencetablefork=2:nfori=k:n%ifabs(X(i+k)-X(i))eps%error('theDATAiserror');%return%endD(i,k)=(D(i,k-1)-D(i-1,k-1))/(X(i)-X(i-k+1));endend%DeterminethecoefficientsoftheNewtoointerpolating%polynomial%计算牛顿插值公式y0=0;fori=1:nz=1;fork=1:i-1z=z*(x0-X(k));endy0=y0+D(i,i)*z;end