function p=ChebyshevInterp(f,n,xx)
n = n-1;                            % Number intervals
Xn = cos((pi/(n+1))*(0.5+(0:n)));   % Chebyshev Nodes
f_X = f(Xn);                        % f at Chebyshev nodes,
omega = BaryWeigths(Xn);            % Barycentric Weigths.
p = zeros(1,length(xx));            % Setting interpolant
is = ismember(xx,Xn);               % Position of a member of Xn at xx.
p(is) = f(xx(is));                  % p(x_k) =f(x_k).
xx = setdiff(xx,Xn);                % xx \setminus Xn.
A = repmat(Xn',1,100);
A = bsxfun(@minus,xx,A);
A = 1./A;
p(~is)=sum(bsxfun(@times,A,(omega.*f_X)'))...
        ./sum(bsxfun(@times,A,omega'));
end

function omega=BaryWeigths(X)
  n=length(X);
  X=X(:);                          % forced to be a column vector
  Omega=repmat(X,1,n);             % [X|X|X|...|X]
  Omega=bsxfun(@minus,X',Omega);   % [x_0-X|x_1-X|...|x_n-X]
  Omega=Omega+eye(n);              % diagonal filled with 1's
  omega=1./prod(Omega);            % Barycentric Weigths
end