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
# on interval [xi, x(i+1)]
#
# 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.');
    return
endif
#
# 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;

endfunction