x=a:(b-a)/n:b; %插值节点
y=f(x);
plot(x,y,'b') %用蓝色线作被插函数图象
hold on
z=a:(b-a)/(2*n):b;
n=length(x);
for j=2:n
for i=n:-1:j
y(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));%计算差商
end
end
u=y(n);
m=length(z);
for j=1:m
for i=n-1:-1:1
u=y(i)+u*(z(j)-x(i)); %计算牛顿插值多项式的值
v(j)=u;
end
u=y(n);
end
plot(z,v,'r') %用红色线作牛顿插值多项式图象
hold off