There are two parts in the statement. First, weak solutions are classical solutions. Second, that they are real analytic.

Weak = classical. This is very simple for linear ODEs, by duality. We know that the classical solutions form an $n$-dimensional space. Then it is enough to prove that the image of ${\mathcal D}(a,b)$ under the adjoint operator $P^*$ is of codimension $n$ at most; this will ensure that the weak solutions of $Pu=0$ form a linear space of dimension $\le n$. Since it contains the space of classical solutions, they are equal.

The codimension $\le n$ is true, because if you take a $\phi\in{\mathcal D}(a,b)$ and solve the Cauchy problem $P^*z=\phi$ with $z(a)=z'(a)=\cdots=z^{n-1}(a)=0$, then $z$ in in ${\mathcal D}(a,b)$ if and only if it vanishes at $b$ together its derivatives up to order $n-1$. Becase these are $n$ linear conditions, we get the codimension $\le n$.

Real analyticity follows from the theory of differential equations with holomorphic coefficients. It has been known since the XIXth century, probably by Cauchy himself. The important point, as in the real case, is that $a_n$ does not vanish. Fuchs went much further by studying regular singularities in holomorphic ODEs, when $a_n$ has simple zeroes. Even the case of irregular singularities (higher order zeroes of $a_n$) is interesting. This yields monodromy groups. It is an important open question whether given points $z_1,\ldots,z_m$ and monodromy groups $G_1,\ldots,G_m$, there exists a complex ODE with singularities $Z_j$ and monodromy group $G_j$ at $z_j$.

There are two parts in the statement. First, weak solutions are classical solutions. Second, that they are real analytic.

Weak = classical. This is very simple for linear ODEs, by duality. We know that the classical solutions form an $n$-dimensional space. Then it is enough to prove that the image of ${\mathcal D}(a,b)$ under the adjoint operator $P^*$ is of codimension $n$ at most; this will ensure that the weak solutions of $Pu=0$ form a linear space of dimension $\le n$. Since it contains the space of classical solutions, they are equal.

The codimension $\le n$ is true, because if you take a $\phi\in{\mathcal D}(a,b)$ and solve the Cauchy problem $P^*z=\phi$ with $z(a)=z'(a)=\cdots=z^{n-1}(a)=0$, then $z$ is in ${\mathcal D}(a,b)$ if and only if it vanishes at $b$ together with its derivatives up to order $n-1$. Because these are $n$ linear conditions, we get the codimension $\le n$.

Real analyticity follows from the theory of differential equations with holomorphic coefficients. It has been known since the XIXth century, probably by Cauchy himself. The important point, as in the real case, is that $a_n$ does not vanish. Fuchs went much further by studying regular singularities in holomorphic ODEs, when $a_n$ has simple zeroes. Even the case of irregular singularities (higher order zeroes of $a_n$) is interesting. This yields monodromy groups, which are subgroups of $GL_n({\mathbb C})$, defined up to conjugacy. It is an important open question whether given points $z_1,\ldots,z_m$ and monodromy groups $G_1,\ldots,G_m$, there exists a complex ODE with singularities $z_j$ and monodromy group $G_j$ at $z_j$.