function [xPlot, pPlot] = cubicSplineNaturalF(t,y,nPlotPts)
% Natural cubic spline interpolation
% Data points are given by (t1,y1), (t2,y2),...
% nPlotPts is the number of points per subinterval to plot.
% Output coefficients are then used to plot the interpolant
% K. Ming Leung, 02/28/03

% t=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3].';
% y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25].';
N=length(t)-1;    % N=20 here
for k=1:N
    h(k)=t(k+1)-t(k);
    b(k)=6*(y(k+1)-y(k))/h(k);
end
u(2)=2*(h(1)+h(2));
v(2)=b(2)-b(1);
for k=3:N
    u(k)=2*(h(k)+h(k-1))-h(k-1)^2/u(k-1);
    v(k)=b(k)-b(k-1)-h(k-1)*v(k-1)/u(k-1);
end
z(N+1)=0;
for k=N:-1:2
    z(k)=(v(k)-h(k)*z(k+1))/u(k);
end
z(1)=0;
for k=1:N
    B(k)=-h(k)*z(k+1)/6-h(k)*z(k)/3+(y(k+1)-y(k))/h(k);
    C(k)=0.5*z(k);
    D(k)=(z(k+1)-z(k))/(6*h(k));
end;

A=y(1:N)';
xPlot=[]; pPlot=[];
for k=1:N
    x=linspace(t(k),t(k+1),nPlotPts);
    p=HornerShift([A(k) B(k) C(k) D(k)],x,t(k));
    xPlot=[xPlot x]; pPlot=[pPlot p];
end
% plot(t,y,'r*',xPlot,pPlot,'b-'); axis equal;
% title('Natural Cubic Spline Interpolation: Dorsal Part of a Duck');
