Dot product with which polynomial gives evaluation at $x_0$?

I thought of following situation/problem, and was surprised that the solution did not jump out at me. Fix a positive integer $n$ and let $V_n$ be the finite-dimensional vector space of all polynomials of degree $\leq n$ with real coefficents. Make $V_n$ into a inner product space with respect to $$ \langle p,q \rangle = \int_{-1}^1 p(x) q(x) \ dx.$$ Now, fix some point $x_0 \in \mathbb{R}$ and consider the linear functional \begin{align*} \varphi_{x_0} : V_n \to \mathbb{R} && \varphi_{x_0}(p) = p(x_0). \end{align*} Because we are in an inner product space, there is a unique polynomial $p_{x_0,n}$ such that \begin{align*} \varphi_{x_0}(p) = \langle p_{x_0,n},p \rangle =\int_{-1}^1 p_{x_0,n}(x) p(x) \ dx && \forall p \in V_n. \end{align*}

Question: What is this polynomial!?

One approach would be to apply the Gramm-Schmidt procedure to $1,x,\ldots,x^n$ to get an orthonormal basis $p_0,p_1,\ldots,p_n$ for $V_n$ and then calculate $\varphi_{x_0}$ on this basis. The required vector should be then $\sum_{i=0}^n \varphi_{x_0}(p_i) p_i$, but I am hoping for a more enlightening representation.

This will be not a very nice and effective answer but gives a more or less direct formula for the polynomial $\varphi_{x_0}$.

We will repeat the proof for Rodriguez's formula for the Legendre polynomials. Let $q_0,q_1,q_2,\ldots,q_{n+1}$ be the sequence of polynomials such that $$ q_0 = \varphi_{x_0}; \qquad q_{k+1}' = q_k \quad\text{and}\quad q_{k+1}(-1)=0 \quad (k=0,1,\ldots,n). $$ We will compute the polynomial $q_{n+1}$.

The condition, applied to the polynomial $(1-x)^k$ (with $0\le k\le n$) yields $$ \int_{-1}^1 \varphi_{x_0}(x) \cdot (1-x)^k dx = (1-x_0)^k. $$ Integrating by parts $k$ times, \begin{align*} (1-x_0)^k &= \int_{-1}^1 q_0(x) \cdot (1-x)^k dx = \underbrace{\Big[q_1(x) \cdot (1-x)^k\Big]_{-1}^1}_0 + \int_{-1}^1 q_1(x) \cdot k(1-x)^{k-1} dx = \\ &= k\int_{-1}^1 q_1(x) \cdot (1-x)^{k-1} dx = k(k-1)\int_{-1}^1 q_2(x) \cdot (1-x)^{k-2} dx = \dots = \\ &= k! \int_{-1}^1 q_k(x) dx = k! \cdot q_{k+1}(1) = k! \cdot q_{n+1}^{(n-k)}(1). \end{align*} Therefore, $$ q_{n+1}^{(n-k)}(1) = \frac{(1-x_0)^k}{k!}. $$

The $n$th Taylor-polynomial of $q_{n+1}$ around $1$ is $$ \sum_{k=0}^n \frac{q_{n+1}^{(n-k)}(1)}{(n-k)!}(x-1)^{n-k} = \sum_{k=0}^n \frac{(1-x_0)^k}{k!}(x-1)^{n-k} = \\ = \frac1{n!} \sum_{k=0}^n \binom{n}{k}(1-x_0)^k(x-1)^{n-k} = \frac{(x-x_0)^n}{n!} $$ so $$ q_{n+1}(x) \equiv \frac{(x-x_0)^n}{n!} \pmod{(x-1)^{n+1}}. \tag1 $$ By the condition $q_1(-1)=q_2(-1)=\ldots=q_{n+1}(-1)=0$ we know that $$ q_{n+1}(x) \equiv 0 \pmod{(x+1)^{n+1}}. \tag2 $$ Since $\deg q_{n+1}\le 2n+1$, by the Chinese Remainder Theorem, the congruences (1) and (2) uniquely determine $q_{n+1}$.

A solution of the congruences (1) and (2), but with higher degree is $$ \left(\bigg(\frac{x-1}2\bigg)^{n+1}-(-1)^{n+1}\right)^{n+1} \cdot \frac{(x-x_0)^n}{n!}. $$ Dividing with remainders, we can see that $$ q_{n+1}(x) = \left(\bigg(\frac{x-1}2\bigg)^{n+1}+(-1)^n\right)^{n+1} \cdot \frac{(x-x_0)^n}{n!} \mod (1-x^2)^{n+1} $$ and therefore $$ \varphi_{x_0}(x) = \left( \left(\bigg(\frac{x-1}2\bigg)^{n+1}+(-1)^n\right)^{n+1} \cdot \frac{(x-x_0)^n}{n!} \mod (1-x^2)^{n+1}\right)^{(n+1)} . $$