function d = natural_spline(x,f)

% function d = natural_spline(x,f)
% derivatives d(k) of natural spline interpolant 
% to data x(k),f(k)

n = length(x);
delta = x(2:n)-x(1:n-1);
df = f(2:n)-f(1:n-1);

A = diag([1./delta(1:n-2); 1/delta(n-1)], -1) + ... 
    diag([2/delta(1); 2./delta(1:n-2) + ... 
          2./delta(2:n-1); 2/delta(n-1)]) + ... 
    diag([1/delta(1); 1./delta(2:n-1)], 1);

b = [3*df./delta.^2; 0] + [0; 3*df./delta.^2];

d = A\b;