The differential operator $D:P_3(R)\rightarrow P_3(R)$ $$D=(x^2-x)\frac{d^2}{dx^2}+(2x-1)\frac{d}{dx}$$ Have a basis of eigenvector, $\Big\{1,x-\frac{1}{2},x^2-x+\frac{1}{6},x^3-\frac{3}{2}x^2+\frac{3}{5}x-\frac{1}{20}\Big\}$
Interestingly, this is also a orthogonal base of the inner product space $P_3(R)$ with the inner product $$\langle p,q\rangle=\int_{0}^{1}p(x)q(x) dx$$
Is there a general relationship between the eigenvector of differential operator of polynomial and orthogonal bases of the inner product space with the above inner product? Why is there such relationship?
The spectral theorem tells you that a self-adjoint linear operator with respect to some inner product has a basis of eigenvectors which is orthogonal in that inner product. Thus it suffices to check that for all $p,q \in P_3$, $\langle p,Dq \rangle=\langle Dp,q \rangle$; then your result follows. To check that you'll need integration by parts; it is crucial here that $(x^2-x)'=2x-1$.
Note that really this is just the finite dimensional spectral theorem from linear algebra, just framed in the language of polynomials, derivatives, and integrals. But it is not really about function spaces. Effectively the same result would happen in $\mathbb{R}^4$ too, if you wrote out $D$ in matrix form.