function [X,Y] = natural_cubic_spline(x,y,N)
% We're given data in the vectors x and y.  The output will be a
% fitting where there are N-1 new points introduced between each
% x(j) and x(j+1). 

% if we have meshpoints x_1...x_n and assume a cubic function on 
% each subinterval [x_j,x_{j+1}] where s''(x_j) =: M_j and s(x_j) = y_j
% then 
%    s(x) = ((x(j+1)-x)^3*M(j) + (x-x(j))^3*M(j+1))/(6*(x(j+1)-x(j)))
%               + ((x(j+1)-x)*y(j) + (x-x(j))*y(j+1))/(x(j+1)-x(j))
%                 -1/6 (x(j+1) - x(j))*((x(j+1) - x)*M(j) + (x-x(j))*M(j+1))
% 
% the requirement that s' be continuous on [x_1,x_n] leads to n-2
% equations for j=2..n-1:
%
% (x_j - x_(j-1))/6 M_{j-1} + (x_{j+1}-x_{j-1})/3 M_j + (x_{j+1}-x_j)/6 M_{j+1}
%                = (y_{j+1}-y_j)/(x_{j+1}-x_j) - (y_j-y_{j-1})/(x_j - x_{j-1})
%
% for the n unknowns M.  We close the system with two more constraints:
% M(1) = 0 and M(n) = 0;

n = length(x);
for j=2:n-1
  A(j,j-1) = (x(j)-x(j-1))/6;
  A(j,j) = (x(j+1)-x(j-1))/3;
  A(j,j+1) = (x(j+1)-x(j))/6;
  b(j) = (y(j+1)-y(j))/(x(j+1)-x(j)) - (y(j)-y(j-1))/(x(j)-x(j-1));
end
% the following will force M(1) = 0 and M(n) = 0 which 
% will create the "natural spline function"
A(1,1) = 1;
b(1) = 0;
A(n,n) = 1;
b(n) = 0;
b = b';
M = A\b;

% now that I have the coefficients M, I can construct the not-a-knot
% spline.  Recall that N is the number of subintervals that will be
% introduced between x(j) and x(j+1).  (There will be N-1 new points.)
j=1;
dx = (x(j+1)-x(j))/N;
xx = x(j):dx:x(j+1)-dx;
X(1:N) = xx;
Y(1:N) = ((x(j+1)-xx).^3*M(j) + (xx-x(j)).^3*M(j+1))/(6*(x(j+1)-x(j))) ...
               + ((x(j+1)-xx)*y(j) + (xx-x(j))*y(j+1))/(x(j+1)-x(j)) ...
                 -(1/6)*(x(j+1) - x(j))*((x(j+1) - xx)*M(j) + (xx-x(j))*M(j+1));   
for j=2:n-1
  dx = (x(j+1)-x(j))/N;
  xx = x(j):dx:x(j+1)-dx;
  X((j-1)*N+1:j*N) = xx;
  Y((j-1)*N+1:j*N) = ((x(j+1)-xx).^3*M(j) + (xx-x(j)).^3*M(j+1))/(6*(x(j+1)-x(j))) ...
               + ((x(j+1)-xx)*y(j) + (xx-x(j))*y(j+1))/(x(j+1)-x(j)) ...
                 -(1/6)*(x(j+1) - x(j))*((x(j+1) - xx)*M(j) + (xx-x(j))*M(j+1));
end
X(N*(n-1)+1) = x(n);
Y(N*(n-1)+1) = y(n);