function [a, b, c, d] = cubic_spline(x,y,nc)
% set nc = 0 for clamped at zero
% = anything else for natural
% Define differences
dy = y(2:end)-y(1:end-1);
h = x(2:end)-x(1:end-1);
% Cubic spline: Natural and clamped use H and u
H = diag(2*(h(1:end-1)+h(2:end)))+diag(h(2:end-1),1)+diag(h(2:end-1),-1);
r = 6./h;
z = zeros(length(r)-1,1);
R = [diag(r(1:end-1)) z z] + [z -diag(r(1:end-1)+r(2:end)) z] + [z z diag(r(2:end))];
u = R*y;
if nc == 0 % Clamped flat (zero 1st deriv at boundary)
 H(1,1) = 1.5*h(1) + 2*h(2);
 H(end,end) = 1.5*h(end) + 2*h(end-1);
 u(1) = u(1) - 3*dy(1)/h(1);
 u(end) = u(end) + 3*dy(end)/h(end);
 mc = H\u;
 m = [3*dy(1)/(h(1)^2)-mc(1)/2;
 mc;
 -3*dy(end)/(h(end)^2)-mc(end)/2];
else % Natural spline (zero 2nd deriv at boundary)
 m = [0; H\u; 0];
end
% To make calculation easier shift m by 1.
me = m(2:end);
m1 = m(1:end-1);

a = y(1:end-1);
b = dy./h - h.*(2*m1 + me)/6;
c = m1/2;
d = (me-m1)./(6*h);
