clear all; close all;

load ducktop.dat
xd=ducktop(:,1);
yd=ducktop(:,2);
plot(xd,yd,'O:')

n=length(xd);
xp=linspace(xd(1),xd(n),100);

% interpolating polynomial
coeff=polyfit(xd,yd,n-1);
yp=polyval(coeff,xp);
plot(xd,yd,'k*',xp,yp,'r')

% our natural cubic spline
yp=spline3natural(xd,yd,xp);
figure(2)
plot(xd,yd,'k*',xp,yp,'r')
hold on

% Matlab's spline command (with default not-a-knot BCs)
yp=spline(xd,yd,xp);
plot(xd,yd,'k*',xp,yp,'b')

% two cubic spline interpolants separated at data point #9:(4.7,2.05)
xp1=linspace(xd(1),xd(9),50);
yp1=spline3natural(xd(1:9),yd(1:9),xp1);
figure(3)
plot(xd(1:9),yd(1:9),'k*',xp1,yp1,'b')
axis([0 14 -2 4])
hold on
xp2=linspace(xd(9),xd(n),50);
yp2=spline3natural(xd(9:n),yd(9:n),xp2);
plot(xd(9:n),yd(9:n),'k*',xp2,yp2,'b')

% quadratic spline
yp=spline22interp(xd,yd,xp);
figure(4)
plot(xd,yd,'k*',xp,yp,'b'), axis([0 14 -2 4])


