function a = pfit(x, y)
	% Compute the interpolating polynomial in Newton form
	% Arguably not the most efficient way of doing this: switching from 
	% an NxN matrix to an N-vector would greatly reduce storage requirements
	n = length(x);
	M = zeros(n,n);
	M(:,1) = y;

	for j=2:n
		for i=1:n-j+1
			M(i,j) = (M(i+1, j-1) - M(i, j-1)) / (x(i+j-1) - x(i));
		end
	end
	a = M(1, :);
end
