[proofplan]
We first check that the Lagrange basis polynomials $\ell_i$ are well-defined and have degree at most $n$. Their defining property is that $\ell_i(x_k)$ equals $1$ when $k=i$ and equals $0$ otherwise. Summing the data values $y_i$ against these basis polynomials gives a polynomial of degree at most $n$ with the required interpolation values. Uniqueness follows because the difference of two such polynomials has degree at most $n$ but has the $n+1$ distinct roots $x_0,\dots,x_n$.
[/proofplan]
[step:Verify that the Lagrange basis polynomials are well-defined and have degree at most $n$]
Fix $i \in \{0,1,\dots,n\}$. Since the nodes $x_0,x_1,\dots,x_n$ are pairwise distinct, for every $j \ne i$ we have $x_i-x_j \ne 0$. Hence each scalar $(x_i-x_j)^{-1}$ is defined in $\mathbb{R}$.
Define the polynomial function
\begin{align*}
\ell_i:\mathbb{R} &\to \mathbb{R}
\end{align*}
by
\begin{align*}
\ell_i(t)=\prod_{\substack{0 \le j \le n, j \ne i}}\frac{t-x_j}{x_i-x_j}.
\end{align*}
If $n=0$, this product has no factors and is therefore the constant polynomial $1$. If $n \ge 1$, it is the product of $n$ linear polynomials multiplied by nonzero real constants. Therefore $\ell_i \in \mathbb{R}[t]$ and $\deg \ell_i \le n$.
[/step]
[step:Compute the values of the basis polynomials at the interpolation nodes]
Let $i,k \in \{0,1,\dots,n\}$. If $k=i$, then
\begin{align*}
\ell_i(x_i)=\prod_{\substack{0 \le j \le n, j \ne i}}\frac{x_i-x_j}{x_i-x_j}=1,
\end{align*}
again using the empty product convention when $n=0$.
If $k \ne i$, then the factor with $j=k$ occurs in the product defining $\ell_i(x_k)$, and that factor is
\begin{align*}
\frac{x_k-x_k}{x_i-x_k}=0.
\end{align*}
Thus $\ell_i(x_k)=0$. Equivalently,
\begin{align*}
\ell_i(x_k)=
\begin{cases}
1, & k=i,
\end{align*}
\begin{align*}
0, & k \ne i.
\end{cases}
\end{align*}
[guided]
The purpose of the polynomials $\ell_i$ is to isolate one interpolation node at a time. We verify this directly. Fix $i,k \in \{0,1,\dots,n\}$. First suppose $k=i$. Substituting $t=x_i$ into the definition gives
\begin{align*}
\ell_i(x_i)=\prod_{\substack{0 \le j \le n, j \ne i}}\frac{x_i-x_j}{x_i-x_j}.
\end{align*}
Every denominator is nonzero because the nodes are pairwise distinct, so each factor equals $1$. Therefore $\ell_i(x_i)=1$. If $n=0$, there are no factors, and the empty product is also $1$.
Now suppose $k \ne i$. The index $j=k$ is included in the product defining $\ell_i$, because that product runs over all $j \ne i$. At $t=x_k$, the corresponding numerator is $x_k-x_k=0$, while the denominator $x_i-x_k$ is nonzero by pairwise distinctness. Thus one factor of the product is
\begin{align*}
\frac{x_k-x_k}{x_i-x_k}=0.
\end{align*}
A product with one zero factor is zero, so $\ell_i(x_k)=0$. Hence $\ell_i$ equals $1$ at its own node and vanishes at all the other interpolation nodes.
[/guided]
[/step]
[step:Construct the interpolating polynomial from the basis polynomials]
Define the polynomial function
\begin{align*}
p:\mathbb{R} &\to \mathbb{R}
\end{align*}
by
\begin{align*}
p(t)=\sum_{i=0}^{n} y_i\ell_i(t).
\end{align*}
Since each $\ell_i$ has degree at most $n$ and each $y_i$ is a scalar, the polynomial $p$ has degree at most $n$.
For each $k \in \{0,1,\dots,n\}$, evaluating at $x_k$ and using the basis values computed above gives
\begin{align*}
p(x_k)=\sum_{i=0}^{n} y_i\ell_i(x_k)=y_k.
\end{align*}
Therefore $p$ is an interpolating polynomial of degree at most $n$.
[/step]
[step:Prove that no other polynomial of degree at most $n$ has the same interpolation values]
We first record the elementary root-count fact needed for uniqueness.
[claim:A nonzero polynomial of degree at most $m$ has at most $m$ distinct real roots]
Let $m \in \mathbb{Z}_{\ge 0}$. If $r \in \mathbb{R}[t]$ is nonzero and $\deg r \le m$, then $r$ has at most $m$ distinct real roots.
[/claim]
[proof]
We use induction on $m$. For $m=0$, a nonzero polynomial of degree at most $0$ is a nonzero constant, so it has no real roots.
Assume the statement holds for $m-1$, where $m \ge 1$, and let $r \in \mathbb{R}[t]$ be nonzero with $\deg r \le m$. If $r$ has no real root, the conclusion holds. Otherwise choose a real root $a \in \mathbb{R}$ of $r$. Write
\begin{align*}
r(t)=\sum_{q=0}^{m} c_q t^q
\end{align*}
with coefficients $c_0,\dots,c_m \in \mathbb{R}$, allowing $c_q=0$ for indices above the actual degree. Since $r(a)=0$, we have
\begin{align*}
r(t)=r(t)-r(a)=\sum_{q=1}^{m} c_q(t^q-a^q).
\end{align*}
For each $q \ge 1$,
\begin{align*}
t^q-a^q=(t-a)\sum_{s=0}^{q-1} t^{q-1-s}a^s.
\end{align*}
Thus there exists $s_0 \in \mathbb{R}[t]$ with $\deg s_0 \le m-1$ such that
\begin{align*}
r(t)=(t-a)s_0(t).
\end{align*}
Because $r$ is nonzero, $s_0$ is nonzero. Every root $b \ne a$ of $r$ satisfies $(b-a)s_0(b)=0$, and since $b-a \ne 0$, it follows that $s_0(b)=0$. By the induction hypothesis, $s_0$ has at most $m-1$ distinct real roots. Therefore $r$ has at most $1+(m-1)=m$ distinct real roots.
[/proof]
Now let $q \in \mathbb{R}[t]$ be another polynomial of degree at most $n$ such that $q(x_i)=y_i$ for every $i \in \{0,1,\dots,n\}$. Define
\begin{align*}
r:\mathbb{R} &\to \mathbb{R}
\end{align*}
by
\begin{align*}
r(t)=p(t)-q(t).
\end{align*}
Then $r \in \mathbb{R}[t]$ and $\deg r \le n$ unless $r$ is the zero polynomial. For every $i \in \{0,1,\dots,n\}$,
\begin{align*}
r(x_i)=p(x_i)-q(x_i)=y_i-y_i=0.
\end{align*}
Thus $r$ has the $n+1$ distinct roots $x_0,x_1,\dots,x_n$. If $r$ were nonzero, the root-count fact with $m=n$ would imply that $r$ has at most $n$ distinct roots, a contradiction. Hence $r$ is the zero polynomial, so $p=q$. This proves uniqueness and completes the proof of the formula.
[/step]