function [b,c,d] = Spcoef(x,f)
%  Calculate coefficients defining a smooth cubic interpolatory spline S.
%
%  Input parameters:
%     x   = vector of n values of the independent variable ordered
%           so that  x(i) < x(i+1)  for all i.
%     f   = vector of values of the dependent variable.
%  Output parameters:
%     b   = vector of S'(x(i)) values.
%     c   = vector of S''(x(i))/2 values.
%     d   = vector of S'''(x(i)+)/6 values (i < n).

b = diff(x);
if any(b <= 0)
  error('Spcoef requires that x(i) < x(i+1) for all i.')
end

n = length(x);

%  Calculate coefficients for the tridiagonal system: store
%  sub-diagonal in b, diagonal in d, difference quotient in c.

c = diff(f) ./ b;
if n == 2
  b(1) = c(1);
  c(1) = 0;
  d(1) = 0;
  b(2) = b(1);
  c(2) = 0;
  return
end
d(1) = 2 * b(1);
for i = 2:n-1
  d(i) = 2 * (b(i) + b(i-1));
end
d(n) = 2 * b(n-1);

%  Calculate estimates for the end slopes.  Use polynomials
%  interpolating data nearest the end.
fp1 = c(1) - b(1) * (c(2) - c(1)) / (b(1) + b(2));
if n > 3
  fp1 = fp1 + b(1) * ((b(1) + b(2)) * (c(3) - c(2)) / ...
       (b(2) + b(3)) - c(2) + c(1)) / (x(4) - x(1));
end      
          
fpn = c(n-1) + b(n-1) * (c(n-1) - c(n-2)) / (b(n-2) + b(n-1));
if n > 3
  fpn = fpn + b(n-1) * (c(n-1) - c(n-2) - (b(n-2) + b(n-1))...
       *(c(n-2) - c(n-3)) / (b(n-2) + b(n-3))) / (x(n) - x(n-3));
end

%  Calculate the right hand side and store it in c.
c(n) = 3 * (fpn - c(n-1));
for k = 2:n-1
  i = n - k + 1;
  c(i) = 3 * (c(i) - c(i-1));
end
c(1) = 3 * (c(1) - fp1);

%  Solve the tridiagonal system.
for k = 2:n
  p = b(k-1) / d(k-1);
  d(k) = d(k) - p * b(k-1);
  c(k) = c(k) - p * c(k-1);
end
c(n) = c(n) / d(n);
for k = n-1:-1:1
  c(k) = (c(k) - b(k) * c(k+1)) / d(k);
end

%  Calculate the coefficients defining the spline.
h = diff(x);
d = diff(c) ./ (3 * h);
b = diff(f) ./ h - h .* (c(1:n-1) + h .* d(1:n-1));
b(n) = b(n-1) + h(n-1) * (2 * c(n-1) + h(n-1) * 3 * d(n-1));
