function output = nspline(x,y,xx)

%NSPLINE  Natural cubic spline interpolation.
%  YY = NSPLINE(X,Y,XX) interpolates the data in vector Y at the 
%  abscissas in vector X using a cubic spline with natural
%  end conditions, and evaluates it at XX.
%  PP = NSPLINE(X,Y) returns the pp-form of the cubic spline 
%  interpolant for later use with PPVAL or PPDER.
%  See also SPLINE, PSPLINE.

% Ref: Kincaid & Cheney, Numerical Analysis, Brooks-Cole, 1991, section 6.4
% Author: Robert Piche, Tampere University of Technology, Finland
%	 (with code from MatLab toolbox routine SPLINE)
% Last updated: Jan. 31, 1994 

x = x(:);
y = y(:);
n = length(x);
[x,ind] = sort(x);
if n<2
  error('There should be at least two data points.')
elseif all(diff(x))==0
  error('The data abscissas should be distinct.')
elseif n~=length(y)
  error('Abscissa and ordinate vectors should be the same length.')
else
  y = y(ind);
  if (n==2), % the interpolant is a straight line
    pp=mkpp(x',[diff(y)./diff(x) y(1)]);
  else
%
% Set up the linear system T*z=v for curvatures (K&C p.319)
%
    h = x(2:n) - x(1:n-1);
    u = 2*( h(1:n-2) + h(2:n-1) );
    b = 6*( y(2:n) - y(1:n-1) )./h;
    if n>3
      T = spdiags([ h(2:n-1) u h(1:n-2)],-1:1,n-2,n-2);
    else
      T = u;
    end
    v = b(2:n-1) - b(1:n-2);
%
% Solve for curvatures at the internal knots
%
    mmdflag = spparms('autommd');   % automatic minimum degree ordering is
    spparms('autommd',0);           % not required for symm. pos. def. T
    z = T\v;
    spparms('autommd',mmdflag);     % reset it to what the user had.
%
% Add the curvatures at the ends
%
    z = [0; z; 0];
%
% Compute coefficients of pp form (K&C Eqn 11)
%
    A = ( z(2:n) - z(1:n-1) )./(6*h);
    B = z(1:n-1)/2;
    C = b/6 - ( z(2:n)/6 + z(1:n-1)/3 ).*h;
    pp = mkpp(x',[ A B C y(1:n-1) ]);
  end
  if nargin==2,
     output = pp;
  else
     output = ppval(pp,xx);
  end
end