function main
% demo1
% demo1b 
 demo2
end

function demo1
 x = [0 1 2 3]
 y = [0 1 4 3]
 clf
 hold on
 splineplot(x,y)
 hold off
end

function demo1b % for comparison with wooden or brass spline
 x = [0 1 2 3]
 y = [0 .3 .0 .1]
 clf
 hold on
 axis equal
 splineplot(x,y)
 hold off
end

function demo2
 x = linspace(-1, 1, 7)
 y = runge(x)
 xx = linspace(-1,1,100);
 yy = runge(xx);
 plot(xx,yy);
 hold on
 splineplot(x,y)
 hold off
end

% Sauer Program 3.6 Cubic spline plot
% Plots spline from data points
% Input: x,y vectors of data points
% Output: none
function splineplot(x,y)
n=length(x);
coeff=splinecoeff(x,y);
            % clear figure window and turn hold on
for i=1:n-1
    x0=linspace(x(i),x(i+1),100);
    dx=x0-x(i);
    y0=coeff(i,3)*dx;   % evaluate using nested multiplication
    y0=(y0+coeff(i,2)).*dx;
    y0=(y0+coeff(i,1)).*dx+y(i);
    plot([x(i) x(i+1)],[y(i) y(i+1)],'o',x0,y0)
end



% Sauer Program 3.5 Calculation of spline coefficients
% Calculates coefficients of cubic spline
% Input: x,y vectors of data points
%   plus two optional extra data v1, vn
% Output: matrix of coefficients b1,c1,d1;b2,c2,d2;...
function coeff=splinecoeff(x,y)
n=length(x);v1=0;vn=0;
A=zeros(n,n);            % matrix A is nxn
r=zeros(n,1);
for i=1:n-1              % define the deltas
    dx(i) = x(i+1)-x(i); dy(i)=y(i+1)-y(i);
end
for i=2:n-1              % load the A matrix
    A(i,i-1:i+1)=[dx(i-1) 2*(dx(i-1)+dx(i)) dx(i)];
    r(i)=3*(dy(i)/dx(i) - dy(i-1)/dx(i-1)); % right-hand side
end
% Set endpoint conditions
A(1,1) = 1;              % natural spline conditions
A(n,n) = 1;
coeff=zeros(n,3);
coeff(:,2)=A\r;          % solve for c coefficients
for i=1:n-1              % solve for b and d
    coeff(i,3)=(coeff(i+1,2)-coeff(i,2))/(3*dx(i));
    coeff(i,1)=dy(i)/dx(i)-dx(i)*(2*coeff(i,2)+coeff(i+1,2))/3;
end
coeff=coeff(1:n-1,1:3);

function r=runge(x)
   r = 1./(1+12*x.^2);