function C = spline_cd ( XI, C )

%*****************************************************************************80
%
%MATLAB version of the FORTRAN SUBROUTINE SPLINE on pages 290/1 of conte/deboor
% With N+1 := length(XI), computes the slope  C(2,i)  at  XI(i) for i=2:N so
% that the piecewise cubic polynomial  sp = CALCCF(Xi,C) has TWO continuous
% derivatives, i.e., is the cubic spline interpolant to the given data.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
[n,m] = size(XI); 
if n>1, XI = reshape(XI,1,n); 
   if m>1, disp(' XI must be a vector'), return, end
end 
[n,m] = size(C);
   if n~=2||m~=length(XI), disp(' XI  and  C  are inconsistent'), return
   end    
N = length(XI)-1;

DX =  diff(XI); 
if any(DX <= 0), disp('entries of XI must be strictly increasing'), return,end
DIAG = [1, diff(C(1,:))./DX]; DX = [0 DX];
C(2,2:N) = 3*(DX(2:N).*DIAG(3:N+1)+DX(3:N+1).*DIAG(2:N));
DIAG(2:N) = 2*(DX(2:N) + DX(3:N+1));
for j=2:N
   G = -DX(j+1)/DIAG(j-1); DIAG(j) = DIAG(j) + G*DX(j-1);
   C(2,j) = C(2,j) + G*C(2,j-1);
end
for j=N:-1:2
   C(2,j) = (C(2,j) - DX(j)*C(2,j+1))/DIAG(j);
end
return, end
% it would be more efficient to compute and output now
% sp = struct('points',XI,'coefs',C);
