% Code 4: Chebyshev Barycentric Interpolation
% Chebyshev nodes of second kind (practical abscissae)
  function ff = barifun(xn,fn,z) 
% Input:  1. xn = interpolation nodes (column vector)  
%         2. fn = function at nodes xn (column vector)
%         3. xn are the Chebyshev nodes mapped to [a,b]:
%            x_j = a + (b-a)*(1+cos(j*pi/n))/2 ; j=0,1,...,n
%         4. z = evaluation grid (z =/= xn ; column vector) 
% Output: ff interpolated function at z
  
   n = length(xn); ff = 0*z; num = ff ; den = ff ; 
   s = (-1).^[0:n-1]';
   sn = s.*[fn(1)/2 ; fn(2:end-1) ; fn(end)/2];
   sd = [s(1)/2 ; s(2:end-1) ; s(end)/2];
   
   for ii = 1:length(num) ;
    num = ((z(ii) - xn).^(-1))'*sn;
    den = ((z(ii) - xn).^(-1))'*sd;
    ff(ii) = num/den;  
   end
