function box=runge(n)
%This routine calculates the table of divided differences
%for the Runge example  y = 1/(x^2 +1), and draws the approx.
xdraw=[-5:0.04:5]';
ydraw=1./(xdraw.^2+1);
plot(xdraw,ydraw,'r')
hold on
%Create the Runge data points.
x=zeros(n,1);
for ii=1:n
    x(ii)=-5+(ii-1)*10/(n-1);
end
y=1./(x.^2+1);
%Calculate the divided difference table.
box=zeros(n,n);
%Initialize the first column of the difference table = y.
for ii=1:n
    box(ii,1)=y(ii);
end
%Compute the table of divided differences.
for k=2:n
    for ii=1:n-k+1
        box(ii,k)=(box(ii+1,k-1)-box(ii,k-1))/(x(ii+k-1)-x(ii));
    end
end
%Now calculate the Newton form of the interpolation polynomial.
term=1+0*xdraw;
poly=box(1,1)*term;
for j=1:n-1
    term=term.*(xdraw-x(j));
    poly=poly+box(1,j+1)*term;
end
plot(xdraw,poly,'b')
hold off