  For $m,n \in \mathbb{N}$, $m\geq 1$ and a given function $f :
\mathbb{R}^m\longrightarrow \mathbb{R}$, the polynomial interpolation problem
(PIP) is to determine a unisolvent node set $P_{m,n} \subseteq \mathbb{R}^m$ of
$N(m,n):=|P_{m,n}|=\binom{m+n}{n}$ points and the uniquely defined polynomial
$Q_{m,n,f}\in \Pi_{m,n}$ in $m$ variables of degree
$\mathrm{deg}(Q_{m,n,f})\leq n \in \mathbb{N}$ that fits $f$ on $P_{m,n}$,
i.e., $Q_{m,n,f}(p) = f(p)$, $\forall\, p \in P_{m,n}$. For $m=1$ the solution
to the PIP is well known. In higher dimensions, however, no closed framework
was available. We here present a generalization of the classic Newton
interpolation from one-dimensional to arbitrary-dimensional spaces. Further we
formulate an algorithm, termed PIP-SOLVER, based on a multivariate divided
difference scheme that computes the solution $Q_{m,n,f}$ in
$\mathcal{O}\big(N(m,n)^2\big)$ time using $\mathcal{O}\big(mN(m,n)\big)$
memory. Further, we introduce unisolvent Newton-Chebyshev nodes and show that
these nodes avoid Runge's phenomenon in the sense that arbitrary periodic
Sobolev functions $f \in H^k(\Omega,\mathbb{R}) \subsetneq
C^0(\Omega,\mathbb{R})$, $\Omega =[-1,1]^m$ of regularity $k >m/2$ can be
uniformly approximated, i.e., $ \lim_{n\rightarrow \infty}||\,f -Q_{m,n,f}
\,||_{C^0(\Omega)}= 0$. Numerical experiments demonstrate the computational
performance and approximation accuracy of the PIP-SOLVER in practice. We expect
the presented results to be relevant for many applications, including numerical
solvers, quadrature, non-linear optimization, polynomial regression, adaptive
sampling, Bayesian inference, and spectral analysis.
