function dvided_diff nodes=[-2,-1,0,1,2]; % nodes n=size(nodes); n=n(2); % determine the number of nodes diff=zeros(n,n); % set up a storage for the divided differences diff(:,1)=(f(nodes))'; % sets up values of the function at the nodes = zeroth difference for j=2:1:n for i=1:1:n-j+1 diff(i,j)=(diff(i+1,j-1)-diff(i,j-1))/(nodes(i+j-1)-nodes(i)); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plot the function versus interpolating polynomial %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1); subplot(1,1,1); clf; % clear the plot % Plot the exact function a=min(nodes); %% figure the left end of the interval b=max(nodes); %% figure the right end of the interval x=[a:.1:b]; plot(x, f(x), '-r'); hold on; % to keep the plot on % calculate the approximated function m=size(x); m=m(2); y = zeros(1,m); for i=1:1:m y(i)=Poly(diff(1,:),nodes,x(i)); end plot(x,y,'-b'); hold off; function P = Poly(coeff, nodes, x) n=size(nodes); n=n(2); % determine the number of nodes %it is expected that size of nodes is the same as the size of coeff! % P=0; for i=1:1:n q=1; for j=1:1:i-1 q=q*(x-nodes(j)); end P=P+coeff(i)*q; end function y=f(x); y=sin(x); %y=3.^x