polyinterp.m

function p = polyinterp(t,y)
% POLYINTERP Polynomial interpolation by the barycentric formula.
% Input:
%   t   interpolation nodes (vector, length n+1)
%   y   interpolation values (vector, length n+1)
% Output:
%   p   polynomial interpolant (function)

t = t(:);                    % column vector
n = length(t)-1;
C = (t(end)-t(1)) / 4;       % scaling factor to ensure stability
tc = t/C;

% Adding one node at a time, compute inverses of the weights.
omega = ones(n+1,1);
for m = 1:n
    d = (tc(1:m) - tc(m+1));    % vector of node differences
    omega(1:m) = omega(1:m).*d; % update previous 
    omega(m+1) = prod( -d );    % compute the new one
end
w = 1./omega;                   % go from inverses to weights

p = @evaluate;

    function f = evaluate(x)
        % % Compute interpolant, one value of x at a time.
        f = zeros(size(x));
        for j = 1:numel(x)
            terms = w ./ (x(j) - t );
            f(j) = sum(y.*terms) / sum(terms);
        end
        
        % Apply L'Hopital's Rule exactly.
        for j = find( isnan(f(:)) )'      % divided by zero here
            [~,idx] = min( abs(x(j)-t) );   % node closest to x(j)
            f(j) = y(idx);                  % value at node
        end        
    end

end