I'm having difficulty figuring out how the bound in the below example is determined.
Consider the function $f(x)=e^{-x}$.
For sample points $\{-1,1\}$, the Lagrange basis interpolants are $\ell_0(x)=-\frac{1}{2}(x-1)$ and $\ell_1(x)=\frac{1}{2}(x+1)$, and the polynomial is $p_1(x)=-e^1\frac{(x-1)}{2}+e^{-1}\frac{(x+1)}{2}$. The error on the interval $[-2,2]$ is bounded by $$|f(x)-p_1(x)|\leq\frac{e^2}{2!}(x-1)(x-2)$$
I understand the bound is supposed to be equal to $$\frac{f^{(N+1)}(\xi(x))}{(N+1)!}\prod_{n=0}^N(x-x_n)$$ but I have no idea why $f^{(N+1)}(\xi(x))=e^2$.
Because the interval for the error estimate is $[-2,2]$. The argument for the estimate is from the full interval, even if the interpolation points $-1,1$ are inside.
See https://proofwiki.org/wiki/Lagrange_Polynomial_Approximation, note that "here exists $ξ$ in the interval spanned by $x$ and $x_i$, $i=0,…,n$" means the convex hull of all $(n+2)$ points. Apart from that, a nicely written proof.