function [a,b,c] = pwQSpline(x,y,muL)
%
% Pre:
%   x,y          column n-vectors. n >= 2 and x(1) < ... < x(n)
%   muL          the value of the first derivative at x(1).
%
% Post:
%   a,b,c        column (n-1)-vectors that define the spline.
%                On [x(i),x(i+1)], the spline S(z) is specified by 
%                the quadratic
%  
%                a(i) + b(i)*(z-x(i)) + c(i)*(z-x(i))^2
%



n = length(x);
a = zeros(n-1,1);
b = zeros(n-1,1);
c = zeros(n-1,1);

% Loop over the intervals left-to-right to compute the interpolant.
% Invariant: derivxi is the value of the first derivative at x(i).

derivxi = muL;
for i = 1 : n - 1
  % If q(x) = a(i) + b(i)* (z-x(i)) + c(i)*(z-x(i))^2
  % then q(x(i)) = a(i), ...

  a(i) = y(i);
  
  % and  q'(x(i))= b(i), ...

  b(i) = derivxi;
  
  % and q(x(i+1)) = a(i)+b(i)*(x(i+1)-x(i))+c(i)*(x(i+1)-x(i))^2,
  
  delta = x(i+1) - x(i);
  c(i) = (y(i+1) - a(i) - b(i) * delta) / delta^2;
  
  % and finally, q'(x(i+1)) = b(i)+2*c(i)*(x(i+1)-x(i))
  % so update derivxi for the next iteration.
  
  derivxi = b(i) + 2*c(i)*delta;
end


