function [b,c,d] = ncspline(x,y)
% Usage: [b,c,d] = ncspline(x,y)
%
% natural cubic spline interpolation of (x,y)
%
% Input:
%	x,y	coordinates to be interpolated
% Output:
%	b,c,d	coefficients of cubic polynomials
% s(x) = yi + bi(x-xi) + ci(x-xi)^2 + di(x-xi)^3
%
% Dependencies:
%	tspdChol.m	tridiagonal spd Cholesky factorization
%	tspdSolve.m	solver for tridiagonal spd systems

n = length(x);
h = x(2:n) - x(1:n-1);			% dx
delta = (y(2:n) - y(1:n-1))./h;		% dy/dx
%
% form the diagonals of the tridiagonal matrix
s = h(2:n-2);	% subdiagonal, symmetric
d = 2*(h(1:n-2) + h(2:n-1));	% main diagonal
%
% decompose the matrix
[d, s, rcond] = tspdChol(d, s);
if (1.0 + rcond == 1.0)
    error('Singularity is detected.');
end
%
% form the right-hand side
r(1:n-2) = delta(2:n-1) - delta(1:n-2);
%
% solve for sigma,
sigma = zeros(size(x));
sigma(2:n-1) = tspdSolve(d, s, r);
% set end points, natural spline
sigma(1) = 0; sigma(n) = 0;
%
% compute the coefficients
b = delta - h.*(sigma(2:n) + 2*sigma(1:n-1));
c = 3*sigma(1:n-1);
d = (sigma(2:n) - sigma(1:n-1))./h;