function yi = qspline(t,y,xi)
% QSPLINE  yi = qspline(t,y,xi)
% Computes not-particularly-natural quadratic spline through points
%              t = t(1),...,t(n+1)
%              y = y(1),...,y(n+1)
% and evaluates at the x-values in xi, to give result yi.   Requires
% t(1) < t(2) < ... < t(n+1) and xi in [t(1),t(n+1)], and finds
% quadratic polynomials in this form:
%     Q_i(x) = y_i + b_i (x - t_i) + c_i (x - t_i)^2
% Try:  >> x=0:3; y=exp(x); xp=0:0.01:3; yq = qspline(x,y,xp);
%       >> plot(x,y,'o',xp,exp(xp),xp,yq)
%       >> legend('data','underlying function','quadratic spline')
% Compare to ncspline.m for natural cubic splines.

t = t(:);  y = y(:);          % force t and y to be column vectors
% check that inputs make sense
if (length(t) ~= length(y)), error('t,y must be vectors of same size'),  end
n = length(t) - 1;
if n <= 0, error('must have at least two points to produce spline'), end
if any(xi < t(1)) | any(xi > t(n+1))
  error('xi must be in interior (interval [t(1),t(n+1)])')
end

h = diff(t);                  % gaps:   h(i) = t(i+1) - t(i)
m = diff(y) ./ h;             % slopes: m(i) = (y(i+1) - y(i)) / h(i)

% directly solve bidiagonal system for coeffs c_i, by forward substitution
c = zeros(size(h));   % c_0 = 0 already
for j=2:n
  c(j) = ( m(j)- m(j-1) - h(j-1)*c(j-1) ) / h(j);
end

% determine other coeffs from c:
b = m - h .* c;

% evaluate quadratics; requires finding interval in which xi lies
xi = xi(:);
yi = zeros(size(xi));
for k = 1:length(xi)
  i = 1;
  while xi(k) > t(i+1)
    i = i + 1;
  end
  % now xi(k) is in [t(i),t(i+1)], so apply Hoerner to evaluate Q_i(x):
  dx = xi(k) - t(i);
  yi(k) = y(i) + dx * (b(i) + dx * c(i));
end
