function [a,b,c,d]=naknotsp(knots,data);
%% function [a,b,c,d]=naknotsp(knots,data);
%% The not-a-knot spline for interpolating data at the knots a=x0<...<x_n=b;
%% Here the knots x_1 and x_n-1 are ``inactive', that means the
%% polynomial pieces joining at those points with all derivatives continuous
%% System of equations is solved using crout3v.m
%%     WARNING: THIS PROGRAM REQUIRES AT LEAST 5 KNOTS
%%  INPUT:
%%    knots - distinct points of interpolation as a row vector
%%    data  - data at the interpolation points as a row vector
%%  OUTPUT:
%%     a - column vector of constant terms for the spline on [x_j,x_j-1] 
%%     b - column vector of coefficients of x for the spline on [x_j,x_j-1] 
%%     c - column vector of coefficients of x^2 for the spline on [x_j,x_j-1] 
%%     d - column vector of coefficients of x^3 for the spline on [x_j,x_j-1] 

% determine lengths of the intervals
n=length(knots)-1;
for j=1:n;
  h(j)=knots(j+1)-knots(j);
end;

% determine the right hand side of the system 
   rhs(2)=(data(1)*h(2)-(h(1)+h(2))*data(2)+h(1)*data(3))/h(2);
   rhs(n)=(data(n+1)*h(n-1)-(h(n-1)+h(n))*data(n)+h(n)*data(n-1))/h(n-1);
for i=3:n-1;
   rhs(i)=3*(data(i+1)*h(i-1)-data(i)*(knots(i+1)-knots(i-1))+data(i-1)*h(i));
   rhs(i)=rhs(i)/(h(i-1)*h(i));
end;

% solve the system using crout3v.m
  dd(1)=2*h(1)*h(2)/3+h(1)^2+h(1)^3/(3*h(2));
  dd(n-1)=2*h(n-1)*h(n)/3+h(n)^2+h(n)^3/(3*h(n-1));
  dd(2:n-2)=2*(knots(4:n)-knots(2:n-2));
  subd(n-2)=(h(n)*h(n-1)-h(n)^3/h(n-1))/3;
  supd(1)=(h(2)*h(1)-h(1)^3/h(2))/3;
  supd(2:n-2)=h(3:n-1);subd(1:n-3)=h(2:n-2);
  c(2:n,1)=crout3v(supd,dd,subd,rhs(2:n));
  
% find the other coefficients
   for j=2:n-1;
     b(j,1)=(data(j+1)-data(j))/h(j) - h(j)*(c(j+1,1)+2*c(j,1))/3;
     d(j,1)=(c(j+1,1)-c(j,1))/(3*h(j));
   end;
a=data(2:n-1)';
b=b(2:n-1,1);
c=c(2:n-1,1);   
d=d(2:n-1,1);   
