function [a,b,c,d] = splinec(x,y)

%    SPLINE  -  Coefficients for cubic spline interpolation.

%       Given data vectors x and y,

%       function [a,b,c,d] = splinec(x,y)

%        creates the coefficients for splines of the form

%      Sk(x) = a(k) + b(k)*(x - x(k)) + c(k)*(x - x(k))^2 + d(k)*(x - x(k))^3


n = max(size(x))-1;
np1 = n+1;

%  First, set up the delta x terms and also obtain the easy coefficients "a"
for k=1:n
  delx(k) = x(k+1)-x(k);

  a(k) = y(k);
end

%  Next, we set up the tridiagonal system to solve for the coefficients "c"

for k=2:n
  km1=k-1;
  DL(k) = delx(km1);
  D(k) = 2*(delx(k)+delx(km1));
  DU(k) = delx(k);
  rhs(k) = 3*(y(k+1) - y(k))/delx(k) - 3*(y(k)-y(k-1))/delx(km1);
end

%  Enter the endpoint conditions here

%  Assuming free spline...
DL(1)=0;
D(1) = 1;
DU(1) = 0;
rhs(1) = 0;

DL(np1)=0;
D(np1) = 1;
DU(np1) = 0;
rhs(np1) = 0;

[DL' D' DU']

%  Now, solve the tridiagonal system

c = tridiag(DL, D, DU, rhs)

%  Finally, go back and get the "b" and "d" values using these "c" values.

c(np1) = 0;

for k=1:n
  d(k) = (c(k+1) - c(k)) / (3*delx(k));
  b(k) = (y(k+1) - y(k))/delx(k) - (c(k+1) + 2*c(k))*delx(k)/3;
end
