function v = natspline(x, y, u)
%NATSPLINE  Natural cubic spline function.
%  v = natspline(x, y, u) finds the piecewise cubic interpolatory
%  spline S(x), with S(x(j)) = y(j), and returns v(k) = S(u(k)).

    n = length(x);
    h = diff(x);               % h(i) = x(i+1) - x(i)
    del = diff(y) ./ h;        % del(i) = (y(i+1) - y(i))/h(i)

    % Diagonals of tridiagonal system

    k = n - 2;
    lowe = zeros(1, k); 
    main = zeros(1, k); 
    uppe = zeros(1, k);

    lowe(1:k-1) = h(2:k);
    main(1:k) = 2*(h(2:k+1) + h(1:k));
    uppe(1:k) = h(1:k);

    % Right-hand side
    r = zeros(1, k);
    r(1:k) = 3*(-del(1:k) + del(2:k+1));

    T = spdiags([lowe' main' uppe'], [-1 0 1], k, k);
    p = T\r';
	c = [0 p' 0];

    % b(i) = (1/h(i))(a(i+1) - a(i)) - (1/3)*h(i)*(2*cj(i) + cj(i+1)); 
	b = del - h.*(2*c(1:end-1) + c(2:end))/3;
    % d(i) = (1/(3*h(i))) * (c(i+1) - c(i)); 
	d =(c(2:end) - c(1:end-1)) ./ (3*h);

    % Find subinterval indices k so that x(k) <= u < x(k+1)

    k = ones(size(u));
    for j = 2:n-1
        k(x(j) <= u) = j;
    end

    % Evaluate spline

    s = u - x(k);
    v = y(k) + s.*(b(k) + s.*(c(k) + s.*d(k)));
end
