function z = ncs(t,y)
% Natural cubic spline interpolation. 
% t : Knots. 
% y : interpolation values. 

n = length(t);
h = diff(t);
b = 6*diff(y)./h;
u = 2*(h(1:end-1)+h(2:end));
v = diff(b);
z = zeros(1,n);

% NB! Notice the indeces: h and b starts with 0, u and v with 1. 

% Gauss elimination. 
for i=2:n-2
    u(i) = u(i)-h(i)^2/u(i-1);
    v(i) = v(i)-h(i)*v(i-1)/u(i-1);
end

% Back solution
z(n)=0;
for i=n-2:-1:1
    z(i+1) = (v(i)-h(i+1)*z(i+2))/u(i);
end
z(1)=0; 



