function [X,Y] = not_a_knot_spline(x,y,z,yz,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).  To close the problem, we provide an additional
% pair of points: z(1) and z(2) where x(1) < z(1) < x(2) and
% x(n-1) < z(2) < x(n).  yz(1) is the value at z(1) and yz(2) is
% the value at z(2).

% 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 a not-a-knot spline, we close the system with two more constraints,
% specified at two points, z_1 in [x_1,x_2] and z_2 in [x_{n-1},x_n]

% test case:
% x = [0 1 2 2.5 3 3.5 4];
% y = x.^3 + 2*x.^2 - x;
% z = [.5 3.75];
% yz =  z.^3 + 2*z.^2 - z;

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 create a not-a-knot spline by requiring correct
% interpolation
xx = z(1);
yy = yz(1);
j = 1;
A(j,j) =  (x(j+1)-xx)^3/(6*(x(j+1)-x(j))) - (1/6)*(x(j+1) - x(j))*((x(j+1) - xx));
A(j,j+1) =  (xx-x(j))^3/(6*(x(j+1)-x(j))) - (1/6)*(x(j+1)-x(j))*(xx-x(j));
b(1) = yy - ((x(j+1)-xx)*y(j) + (xx-x(j))*y(j+1))/(x(j+1)-x(j));
xx = z(2);
yy = yz(2);
j = n-1;
A(j+1,j) =  (x(j+1)-xx)^3/(6*(x(j+1)-x(j))) - (1/6)*(x(j+1) - x(j))*((x(j+1) - xx));
A(j+1,j+1) =  (xx-x(j))^3/(6*(x(j+1)-x(j))) - (1/6)*(x(j+1)-x(j))*(xx-x(j));
b(j+1) = yy - ((x(j+1)-xx)*y(j) + (xx-x(j))*y(j+1))/(x(j+1)-x(j));
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);