function a = polyinterp(x, y)
% Returns the length(x) coefficients of a polynomial of degree length(x)-1
% interpolating the points (x,y).

  n = length(x);
  if length(y) ~= n, stop('x and y must be of equal length'), end
  V = zeros(n, n);
  for k = 1:n
      V(:,k) = x .^ (n-k);
  end
  a = V \ y;    % try different methods for solving the system
  
  return