% chebeval(f, xx)
%
% Evaluate a Chebyshev interpolant by barycentric interpolation.
% See Trefethen's "Approximation Theory and Approximation Practice".
% http://www2.maths.ox.ac.uk/chebfun/ATAP/

function fxx = chebeval(f, xx)

N = length(f)-1;
fxx = 0*xx;

for i = 1:length(xx)
done = 0;

z = xx(i);
zd = (z-1);
if zd == 0
fxx(i) = f(1);
done = 1;
continue;
end
pnum = f(1)/zd/2;
pden = 1/zd/2;

for j = 1:N-1
xj = cos(pi*j/N);
zd = (z-xj);
if zd == 0
fxx(i) = f(j+1);
done = 1;
break;
end
pnum = pnum + (-1)^j*f(j+1)/zd;
pden = pden + (-1)^j/zd;
end

zd = (z+1);
if zd == 0
fxx(i) = f(end);
done = 1;
continue;
end
pnum = pnum + (-1)^N*f(end)/zd/2;
pden = pden + (-1)^N/zd/2;

if ~done
fxx(i) = pnum/pden;
end
end