Let $P = a_n(x) D_x^n + a_{n-1}(x) D_x^{n-1} + \ldots + a_0(x)$ be a linear ordinary differential operator with polynomial (or real analytic) coefficients $a_j(x)$. Suppose that $a_n(x)$ doesn't vanish on an interval $(a,b)$ and that $u$ is a weak solution of $P u = 0$ on $(a,b)$. It can be concluded from the elliptic regularity theorem that $u$ is in fact real analytic and a classical solution to the differential equation.

In other words, the operator $P$ is analytically hypoelliptic.

However, appealing to the elliptic regularity theorem seems a bit much for the above case of an ordinary differential operator. Does anyone know a (canonical) reference for this (probably much more classical) case? I browsed through several books on differential equations but the closest I could find was in G. Folland's Fourier analysis and its applications where he mentions this fact (with just $C^\infty$ smoothness) on the top of page 344 without reference or proof.

Update: Thanks for the answers so far! While very helpful, I think I didn't make it clear what I was looking for. For an audience which may not be familiar with the theory of elliptic operators, what would be the proper way to (in one sentence plus a reference) justify that a weak solution in the case at hand is automatically real analytic? Currently, I have to make reference to a textbook on PDEs which treats elliptic regularity even though the problem for linear ODEs seems so much simpler.

Thank you!