function moments=RungeNatSpline(n)
%This routine calculates the second-derivative moments
%for the Runge example  y = 1/(x^2 +1), with equally-spaced
%knot points, and draws the natural-spline approx.
xdraw=[-5:0.04:5]';
ydraw=1./(xdraw.^2+1);
plot(xdraw,ydraw,'r')
hold on
%Create the Runge data points.
hstep=10/(n-1);
x=zeros(n,1);
for ii=1:n
    x(ii)=-5+(ii-1)*hstep;
end
y=1./(x.^2+1);
%Calculate the natural-spline moments.  We ought to do a proper
%tri-diagonal solver, but I can't be bothered.
moments=zeros(n,1);
A=zeros(n-2,n-2);
%Give the natural-spline conditions on the end moments.
moments(1)=0;
moments(n)=0;
%Solve the spline moment conditions.
for ii=1:n-2
    if ii>1
        A(ii-1,ii)=1;
    end
    A(ii,ii)=4;
    if ii<n-2
        A(ii+1,ii)=1;
    end
end
b=zeros(n-2,1);
for ii=2:n-1
    b(ii-1)=6*(y(ii+1)-2*y(ii)+y(ii-1))/(hstep^2);
end
tempM=A\b;
for ii=1:n-2
    moments(ii+1)=tempM(ii);
end
%Now calculate the cubic-spline approximation to the function.
nsize=size(xdraw);
ndraw=nsize(1);
splin=zeros(ndraw,1);
for iplot=1:ndraw
    for ii=1:n-1
        if xdraw(iplot)>=x(ii)
            if xdraw(iplot)<=x(ii+1)
                jval=ii;
            end
        end
    end
    x1min=x(jval+1)-xdraw(iplot);
    xmin=xdraw(iplot)-x(jval);
    term1=(moments(jval)*x1min^3+moments(jval+1)*xmin^3)/(6*hstep);
    term2=(y(jval)/hstep-moments(jval)*hstep/6)*x1min;
    term3=(y(jval+1)/hstep-moments(jval+1)*hstep/6)*xmin;
    splin(iplot)=term1+term2+term3;
end
plot(xdraw,splin,'b')
hold off