数值分析实验，站在巨人（百度）的肩膀上写了代码~~~~~。

function [p,q]=chashang_newton(x,y)
%%%%%%%%%%%%%%%%%%%%差商表的求法 
 n=length(x);
 p(:,1)=x;
 p(:,2)=y;
 for j=3:n+1
    p(1:n+2-j,j)=diff(p(1:n+3-j,j-1)) ./(x(j-1:n)-x(1:n+2-j));
 end
 disp(p);
 q=p(1,2:n+1)';

%%%%%%%%%%%%%%%%%%%%牛顿四阶插值多项式的求法和得到相应的函数值
x1=input('输入节点坐标x=')
y=input('输入节点坐标函数值f(x)=')
x2=input('输入所要计算的节点x2=')
syms x
n=length(x1);
for i=2:n
   f1(i,1)=(y(i)-y(i-1))/(x1(i)-x1(i-1));
end


for i=2:n
   for j=i+1:n
       f1(j,i)=(f1(j,i-1)-f1(j-1,i-1))/(x1(j)-x1(j-i));
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%Newton插值函数
f1=[y',f1]
Newton=f1(1,1);
for i=2:n
   tt=1;
   for j=1:i-1
       tt=tt*(x-x1(j));
   end
   Newton=Newton+f1(i,i)*tt;
end
fprintf('Newton插值函数为n')
expand(Newton)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=x2;
p=eval(Newton);
fprintf('Newton插值函数在所求点x2的函数值为n')
p









%%%%%%%%%%%%%%%%%%%%%%%%%%四阶拉格朗日插值多项式
%输出的量：n次拉格朗日插值多项式L及其系数向量C，基函数l及其系数矩阵L
%本程序为Lagrange1插值，其中x1,y1%为插值节点和节点上的函数值，输出为插值点xx的函数值，%xx可以是向量。
x2=input('本程序为Lagrange1插值，其中x2插值节点:')
y2=input('本程序为Lagrange1插值，其中y2节点上的函数值:')
xx=input('本程序为Lagrange1插值，其中x1插值节点:')
syms x
n=length(x2);
for i=1:n
t=x2;t(i)=[];L(i)=prod((x-t)./(x2(i)-t));% L向量用来存放插值基函数
end
u=sum(L.*y2);
p=simplify(u) % p是简化后的Lagrange插值函数（字符串）
yy=subs(p,x,xx);    %p是以x为自变量的函数，并求xx处的函数值

%%%%%%%%%%%%%%%%%%%%牛顿插值多项式和拉格朗日插值多项式图像
t=input('请输入更多插值点例如：t=[0.5,2.3,3.9]：-------------')
x=input('牛顿插值多项式和拉格朗日插值多项式图像，其中x插值节点:')
y=input('牛顿插值多项式和拉格朗日插值多项式图像，其中y节点上的函数值:')
syms p;
n=length(x);
s=0;
for(k=1:n)
    la=y(k);
    for(j=1:k-1)
        la=la*(p-x(j))/(x(k)-x(j));
    end;
    for(j=k+1:n)
        la=la*(p-x(j))/(x(k)-x(j));
    end;
    s=s+la;
    simplify(s);
end
if nargin==2,
        s=subs(s,'p','x');
        s=collect(s);
        s=vpa(s,6);
    else
        m=length(t);
        for i=1:m,
            temp(i)=subs(s,'p',t(i));
        end
        s=temp;
        yh=lag(x,y)
        ezplot(yh,[-4,5])
end
yh=lag(x,y)
ezplot(yh,[-4,5])

function s=lag(x,y,t)
syms p;
n=length(x);
s=0;
for(k=1:n)
    la=y(k);
    for(j=1:k-1)
        la=la*(p-x(j))/(x(k)-x(j));
    end;
    for(j=k+1:n)
        la=la*(p-x(j))/(x(k)-x(j));
    end;
    s=s+la;
    simplify(s);
end
if nargin==2,
        s=subs(s,'p','x');
        s=collect(s);
        s=vpa(s,6);
    else
        m=length(t);
        for i=1:m,
            temp(i)=subs(s,'p',t(i));
        end
        s=temp;
        yh=lag(x,y)
        ezplot(yh,[-4,5])
end
    

 x=[1,2,3,-4,5]'

x =

     1
     2
     3
    -4
     5

>> y=[2,48,272,1182,2262]'

y =

           2
          48
         272
        1182
        2262

>> [p,q]=chashang_newton(x,y)
           1           2          46          89           6           4
           2          48         224          59          22           0
           3         272        -130         125           0           0
          -4        1182         120           0           0           0
           5        2262           0           0           0           0

输入节点坐标x=[1,2,3,-4,5]

x1 =

     1     2     3    -4     5

输入节点坐标函数值f(x)=[2,48,272,1182,2262]

y =

           2          48         272        1182        2262

输入所要计算的节点x2=-4

x2 =

    -4


f1 =

           2           0           0           0           0
          48          46           0           0           0
         272         224          89           0           0
        1182        -130          59           6           0
        2262         120         125          22           4

Newton插值函数为n 
ans =
 
4*x^4 - 2*x^3 + x^2 - 3*x + 2
 
Newton插值函数在所求点x2的函数值为n
p =

        1182

本程序为Lagrange1插值，其中x2插值节点:[1,2,3,-4,5]

x2 =

     1     2     3    -4     5

本程序为Lagrange1插值，其中y2节点上的函数值:[2,48,272,1182,2262]

y2 =

           2          48         272        1182        2262

本程序为Lagrange1插值，其中x1插值节点:-4

xx =

    -4

 
p =
 
4*x^4 - 2*x^3 + x^2 - 3*x + 2