function S = ncseval(t,y,z,x)
n = length(t);
S = zeros(1,length(x));
h = diff(t);

for i=1:n-1
    % Find the x's in (t(i),t(i+1)].
    ni = find(x>=t(i) & x < t(i+1));
    if i==n-1 
        ni = [ni,find(x==t(end))];
    end
    % Calculate S.
    S(ni) = (z(i)*(t(i+1)-x(ni)).^3+z(i+1)*(x(ni)-t(i)).^3)/(6*h(i)) ...
          + (y(i+1)/h(i)-z(i+1)*h(i)/6)*(x(ni)-t(i)) ...
          + (y(i)/h(i)-z(i)*h(i)/6)*(t(i+1)-x(ni));
end
