function yp = polyint(x,y,xp)
% function polyint(x,y,xp)
%   perform polynomial interpolation
%   input:  x,y - data vectors
%           xp  - new abscissa vector 
%   output: yp  - polynomial interpolation on xp value(s)
%

n = length(x);
for i=1:n
  for j=1:n
    a(i,j) = x(i)^(j-1);
  end
end
b=a\y';
b = b(n:-1:1);      % reverses order of polynomial coefficients
yp = polyval(b,xp); % evaluate polynomial interpolation on xp