function[d]=divdif(x,f) % % Computes divided difference table. x is the vector of interpolation % points, f the vector of data. % n=length(x)-1; d=f; for i=2:n+1 for j=n+1:-1:i d(j)=(d(j)-d(j-1))/(x(j)-x(j-i+1)); end end function p=interp(x,d,t) % % Computes interpolating polynomial at points in the vector t % where x is the vector of nodes and d is the vector of divided % differences produced by divdif. % n=length(x)-1; k=length(t); p=d(n+1)*ones(1,k); for i=n:-1:1 p=d(i)+(t-x(i)).*p; end >> x=[0 1 2 3 4] x = 0 1 2 3 4 >> y=[2 4 8 10 6] y = 2 4 8 10 6 >> d= divdif(x,y) d = 2.0000 2.0000 1.0000 -0.6667 0 >> z=[1 1.5 2 2.5 3 3.5 4] z = 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 >> p=interp(x,d,z) p = 4.0000 6.0000 8.0000 9.5000 10.0000 9.0000 6.0000