function ret = quadraticspline(x,y,xx)
%
% ret = quadraticspline(x,y)
% ret = quadraticspline(x,y,xx)
%
% if (x(1) == x(2)) then y(1) is the value of the derivative at x(1)
% if (x(end-1) == x(end)) then y(end) is the value of the derivative at x(end)
%
x = x(:);
y = y(:);
n = length(x)-1;
if (x(1) == x(2))
% left derivative  
  h = diff(x(2:n+1));
  yder = y(1);
  y = y(2:n+1);
  x = x(2:n+1);
  A = spdiags(ones(n,2),[-1,0],n,n);
  b = [yder;2*diff(y)./h];
else % x(n) == x(n+1)
  h = diff(x(1:n));
  yder = y(n+1);
  y = y(1:n);
  x = x(1:n);
  A = spdiags(ones(n,2),[0,1],n,n);
  b = [2*diff(y)./h;yder];
end
m = A\b;
coefs = [diff(m)./(2*h),m(1:n-1),y(1:n-1)];
ret = mkpp(x,coefs);
if (nargin == 3)
  ret = ppval(ret,xx);
end