function [yinterp, coeffs] = cubicSpline(x, y, u, type, s1, sn)
%
%function [yinterp, coeffs] = cubicSpline(x, y, u, type, s1, sn)
%
% x, y are vectors containing abcissas, ordinates of knot points
% u is vector of abcissas at which to interpolate/extrapolate
% type may be
%   1 for "natural"    (default)
%   2 for "specSlope"  (should be tested)
%   3 for "parabolic"  (should be tested)
%   4 for "notKnot"
% s1, sn
%   These should only be provided if type 2 is used, in which case
%   the values of the cubic spline's 2nd derivative at the first
%   and last data points will be set to s1 and sn respectively.
%
% return values
%   yinterp: a vector of ordinates corresponding to u
%   coeffs:  a matrix, the rows of which consist of the 4 coeffs.
%            a, b, c, d of a cubic polynomial
%               a(x - x_i)^3 + b(x - x_i)^2 + c(x - x_i) + d
%            that fits the two data points for subinterval i.

  if (nargin < 4)
%   type = 'natural';
    type = 1;
  end

  n = length(x);
  A = zeros(n-2, n-2);
  dx = diff(x);
  dy = diff(y);

% if (type == 'notKnot')
  if type == 4
    yinterp = spline(x, y, u);
	return
% elseif (type == 'parabola')
  elseif type == 3
	if n == 3
	  s = 6*diff(dy./dx)/(3*(dx(1) + dx(2)));
	  s = [s; s; s];
	else
      for i = 1:(n-2)
	    if i == 1
	      A(1,1) = 3*dx(1) + 2*dx(2);
	      A(1,2) = dx(2);
	    elseif i == n - 2
	      A(n-2,n-3) = dx(n-2);
	      A(n-2,n-2) = 3*dx(n-1) + 2*dx(n-2);
	    else
		  A(i,i-1) = dx(i);
		  A(i,i) = 2*(dx(i)+dx(i+1));
		  A(i,i+1) = dx(i+1);
	    end
	  end
	  s = A\(6*diff(dy./dx));
	  s = [s(1); s; s(n-2)];
	end
% elseif (type == 'specSlope')
  elseif type == 2
	if n == 3
	  s = (6*diff(dy./dx) - dx(1)*s1 - dx(2)*sn)/(2*(dx(1) + dx(2)));
	  s = [s1; s; sn];
	else
      for i = 1:(n-2)
	    if i == 1
	      A(1,1) = 2*(dx(1) + dx(2));
	      A(1,2) = dx(2);
	    elseif i == n - 2
	      A(n-2,n-3) = dx(n-2);
	      A(n-2,n-2) = 2*(dx(n-1) + dx(n-2));
	    else
		  A(i,i-1) = dx(i);
		  A(i,i) = 2*(dx(i)+dx(i+1));
		  A(i,i+1) = dx(i+1);
	    end
	  end
	  rhs = 6*diff(dy./dx);
	  rhs(1) = rhs(1) - dx(1)*s1;
	  rhs(n-1) = rhs(n-1) - dx(n-1)*sn;
	  s = A\rhs;
	  s = [s1; s; sn];
	end
  else
	if n == 3
	  s = [0; 6*diff(dy./dx)/(2*(dx(1) + dx(2))); 0];
	else
      for i = 1:(n-2)
	    if i == 1
	      A(1,1) = 2*(dx(1) + dx(2));
	      A(1,2) = dx(2);
	    elseif i == n - 2
	      A(n-2,n-3) = dx(n-2);
	      A(n-2,n-2) = 2*(dx(n-1) + dx(n-2));
	    else
		  A(i,i-1) = dx(i);
		  A(i,i) = 2*(dx(i)+dx(i+1));
		  A(i,i+1) = dx(i+1);
	    end
	  end
	  s = [0; A\(6*diff(dy./dx)); 0];
	end
  end

  a = (1/6)*diff(s)./dx;
  b = s(1:n-1)/2;
  c = dy./dx - dx.*(2*s(1:n-1) + s(2:n))/6;
  
% Find subinterval indices k so that x(k) <= u < x(k+1)

  k = ones(size(u));
  for j = 2:n-1
     k(x(j) <= u) = j;
  end

% Evaluate spline

  s = u - x(k);
  yinterp = y(k) + s.*(c(k) + s.*(b(k) + s.*a(k)));
% yinterp = [a'; b'; c'; y(1:n-1)'];
  coeffs = [a'; b'; c'; y(1:n-1)']';