function p = chebint(fk, x)

%  The function p = chebint(fk, x) computes the polynomial interpolant
%  of the data (xk, fk), where xk are the Chebyshev nodes.  
%  Two or more data points are assumed.
%
%  Input:
%  fk:  Vector of y-coordinates of data, at Chebyshev points 
%       x(k) = cos((k-1)*pi/(N-1)), k = 1...N.
%  x:   Vector of x-values where polynomial interpolant is to be evaluated.
%
%  Output:
%  p:    Vector of interpolated values.
%
%  The code implements the barycentric formula; see page 252 in
%  P. Henrici, Essentials of Numerical Analysis, Wiley, 1982.
%  (Note that if some fk > 1/eps, with eps the machine epsilon,
%  the value of eps in the code may have to be reduced.)

%  J.A.C. Weideman, S.C. Reddy 1998

  fk = fk(:); x = x(:);                    % Make sure data are column vectors.

   N = length(fk); 
   M = length(x);
     
  xk = sin(pi*[N-1:-2:1-N]'/(2*(N-1)));    % Compute Chebyshev points.

   w = ones(N,1).*(-1).^[0:N-1]';          % w = weights for Chebyshev formula
w(1) = w(1)/2; w(N) = w(N)/2;
 
   D = x(:,ones(1,N)) - xk(:,ones(1,M))';  % Compute quantities x-x(k)
   D = 1./(D+eps*(D==0));                  % and their reciprocals.
  
   p = D*(w.*fk)./(D*w);                   % Evaluate interpolant as
                                           % matrix-vector products.