I've successfully shown that Legendre polynomials are an orthonormal basis for functions. However, I'm wondering how to proof that all functions f(x) can be written in the form (Gauss-Hermite series):
where $a$ are constants are $e$ are the Legendre polynomials.

Firstly note that $\psi_n (x) = \frac{1}{\sqrt{2^n n!}} \frac{1}{\sqrt[4]{\pi}} e^{-\frac12 x^2}e_n(x)$, $n \in \mathbb N$ are all the eigenfunctions of the quantum harmonic oscillator: $$ - \frac12 \psi_n ''(x) +\frac12 x^2 \psi_n(x) = E_n\psi_n(x), \qquad E_n=n+\frac12. $$ We know that such functions constitute a basis for $L^2(\mathbb R)$, and your result follows for $f \in L^2(\mathbb R)$.