x=x(:); y=y(:); xi=xi(:); n=length(x);
if n~=length(y), error('input x and y arrays must be of same length'), end
dx=diff(x);
if any(dx <= 0)
error('input x-array must be in strictly ascending order');
end
if any(xi<x(1)) || any(xi > x(n))
warning('All interpolation points xi must lie between x(1) and x(n)')
end
m=diff(y)./dx;
mm=2*m(1)-m(2); mmm=2*mm-m(1); % augment at left
mp=2*m(n-1)-m(n-2); mpp=2*mp-m(n-1); % augment at right
m1=[mmm; mm; m; mp; mpp]; % slopes
dm=abs(diff(m1)); f1=dm(3:n+2); f2=dm(1:n); f12=f1+f2;
id=find(f12 > 1e-8*max(f12)); b=m1(2:n+1);
b(id)=(f1(id).*m1(id+1)+f2(id).*m1(id+2))./f12(id);
c=(3*m-2*b(1:n-1)-b(2:n))./dx;
d=(b(1:n-1)+b(2:n)-2*m)./dx.^2;
% Loop replaced by vector ops following tip from johannes.korsawe@volkswagen.de
% 1/19/2006
%
[ncnt,bin]=histc(xi,x);
bin=min(bin,n-1);
bb=bin(1:length(xi));
wj=xi-x(bb);
yi=((wj.*d(bb) +c(bb)).*wj+b(bb)).*wj+y(bb);