function S=cptsfit(X,Y)

% Cubic Parabolically Terminated Spline
%Input    - X is the 1xn abscissa vector
%         - Y is the 1xn ordinate vector

%Output   - S: rows of S are the coefficients from high order to low one for the cubic parabolically terminated spline interpolants

%  NUMERICAL METHODS: Matlab Programs
% (c) 2007 by Xie Liling
%  Complementary Software to accompany the textbook:
%  Information and Computing Science: A Laboratory Course

N=length(X)-1;
H=diff(X);
D=diff(Y)./H;
A=H(2:N-1);
B=2*(H(1:N-1)+H(2:N));
C=H(2:N);
C=H(2:N);
U=6*diff(D);

% Parabolically terminated spline endpoint constraints

B(1)=B(1)+H(1);
B(N-1)=B(N-1)+H(N);

for k=2:N-1
   temp=A(k-1)/B(k-1);
   B(k)=B(k)-temp*C(k-1);
   U(k)=U(k)-temp*U(k-1);
end

M(N)=U(N-1)/B(N-1);

for k=N-2:-1:1
   M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
end

% Parabolically terminated spline endpoint constraints

M(1)=M(2);
M(N+1)=M(N);

for k=0:N-1
   S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));
   S(k+1,2)=M(k+1)/2;
   S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
   S(k+1,4)=Y(k+1);
end