the last matrix before the appendix is not the identity matrix, right now the matrix is: \begin{bmatrix} 1 & 0 & 0 & \dots & 0\\ 1 & 0 & 0 & \dots & 0\\ 1 & 0 & 0 & \dots & 0\\ \vdots & \vdots & \vdots & \ddots &\vdots \\ 1 & 0 & 0 & \dots & 1 \end{bmatrix}

The Lagrange polynomials form the normal basis of most Finite Elements Method (FEM) software. There are other polynomials which are used as well, but these are the workhorse of most solvers.

In the Polynomial Interpolation Theorem, you have r(x) = p(x) - r(x), but I think it should be q(x) = p(x) - r(x).