function [a,b,c,d]=natsp3_4(knots,data);
%% function [a,b,c,d]=natsp3_4(knots,data);
%% The natural spline for interpolating data at the knots a=x0<...<x_n=b;
%% System of equations is solved using crout3v.m
%%    WARNING: THIS PROGRAM REQUIRES AT LEAST 4 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 (without first and last eqn)
for i=2:n;
   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
  c(1,1)=0;
  c(n+1,1)=0;
  dd=2*(knots(3:n+1)-knots(1:n-1));subd=h(1:n-2);supd=h(2:n-1);
  c(2:n,1)=crout3v(supd,dd,subd,rhs(2:n));
  
% find the other coefficients
   for j=1:n;
     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(1:n)';
c=c(1:n,1)   
