function ret = spline_notaknot (x, y, xx)
%
% function pp = spline_notaknot (x, y)
% function yy = spline_notaknot (x, y, xx)
%
% Simula il comportamento di spline.
x = x(:);
[x, idx] = sort (x);
y = y(:);
y = y(idx);
n = length (x);
h = diff (x);
d0 = [h(2:n - 3) / 6; (h(n-2)-h(n-1))/6; 0];
d1 = [(h(1)+2*h(2))/6; (h(2:n-3)+h(3:n-2)) / 3; (2*h(n-2)+h(n-1))/6];
d2 = [0; (h(2)-h(1))/6; h(3:n-2)/6];
A = spdiags ([d0, d1, d2], -1:1, n - 2, n - 2);
d = (y(3:n)-y(2:n-1))./h(2:n-1)-(y(2:n-1)-y(1:n-2))./h(1:n-2);
d(1) = h(2)/(h(1)+h(2)) * d(1);
d(n - 2) = h(n-2)/(h(n-2)+h(n-1)) * d(n - 2);
m = zeros (n, 1);
m(2:n-1) = A \ d;
m(1) = m(2) + h(1) / h(2) * (m(2) - m(3));
m(n) = m(n-1) + h(n-1) / h(n - 2) * (m(n-1) - m(n-2));
coefs(:,1) = (m(2:n)-m(1:n-1))./(6*h(1:n-1));
coefs(:,2) = m(1:n-1)/2;
coefs(:,3) = (y(2:n)-y(1:n-1))./h(1:n-1)-...
             m(2:n).*h(1:n-1)/6-m(1:n-1).*h(1:n-1)/3;
coefs(:,4) = y(1:n-1);
ret = mkpp (x, coefs);
if (nargin == 3)
   ret = ppval (ret, xx);
end
%!test
%! x = rand (1, 4);
%! y = rand (1, 4);
%! assert (spline (x, y), spline_notaknot (x, y), -1000 * eps)
%!test
%! x = rand (1, 5);
%! y = rand (1, 5);
%! assert (spline (x, y), spline_notaknot (x, y), -1000 * eps)
