function v = natsplfun(x, y, b, c, d, u)
%NATSPLFUN  Natural cubic spline function.
%  v = natsplfun(x, y, b, c, d, u) returns v(k) = S(u(k)).

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

    n = length(x); 
    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
