I'm teaching an honors differential equations class and have been using linear algebra heavily. I thought it would be interesting to include a proof of the method of undetermined coefficients along the lines described below, which I figured is already known. Nonetheless, Google turns up no proof other than the obvious one by integration, which is significantly messier and uses no linear algebra, so perhaps I'm wrong.

What is the (standard, best, earliest) reference to the following proof?

Note: I have made a conscious decision not to ask this on math.stackexchange since, based on my last question asking about a reference there, I don't think it will even be noticed.

**Theorem:** Let $L$ be a constant-coefficient $n$th order linear differential operator with characteristic polynomial $p(\lambda)$, and let $y(t) = t^m e^{qt}$ for some integer $m \geq 0$ and $q \in \mathbb{C}$. Suppose that $q$ is a $k$-fold root of $p(\lambda)$. Then the equation $Lx = y$ has a solution of the form $t^k f(t) e^{qt}$, where $f(t)$ is a polynomial of degree $m$.

**Proof outline:**

$y$ is a solution of the equation $L_2 y = 0$, where $L_2 = (\tfrac{d}{dt} - q)^m$. The pair of equations $Lx = y, L_2 y = 0$ can be written as a system of $n + m + 1$ first-order equations via the linear transformation $\vec{z} = (x, x', \dots, x^{(n - 1)}, y, y', \dots, y^{(m)})$, for which the coefficient matrix is: $$M = \left(\begin{array}{c|c} A & \begin{smallmatrix} 0 & \cdots & 0 \\\ \vdots & & \vdots \\\ 1 & \cdots & 0 \end{smallmatrix} \\\ \hline 0 & C \end{array}\right)$$ where $A$ is the companion matrix of $L$ and $C$ is the companion matrix of $L_2$.

It is easily verified that:

The eigenvalues of $M$ are the same as those of $A$, with an additional $m + 1$ repetitions of $q$; each eigenvalue $\lambda$ has a one-dimensional eigenspace, spanned by $(1, \lambda, \dots, \lambda^{n - 1}, p(\lambda), \lambda p(\lambda), \dots, \lambda^m p(\lambda))$.

Every Jordan chain for $A$, namely, vectors $\vec{v}_i$ associated with an eigenvalue $\lambda$ such that $(A - \lambda I) \vec{v}_i = \vec{v}_{i - 1}$, upgrades to vectors $\vec{u}_i$ with the same property for $M$, by attaching $m + 1$ zeros to the end.

It follows that $M$ has a generalized eigenbasis consisting of upgrades of a generalized eigenbasis for $A$, together with $m + 1$ additional generalized eigenvectors filling out a Jordan chain of length $k + m + 1$ for eigenvalue $q$.

A generalized eigenvector $\vec{v}_h$ of eigenvalue $\lambda$ at the top of a Jordan chain produces a fundamental solution to $\vec{z}' = M\vec{z}$ via the matrix exponential: $$e^{\lambda t} \sum_{i = 0}^{h - 1} \frac{1}{i!} t^i \vec{v}_{h - i}.$$ The general solution is therefore the general linear combination of these. Taking the first components (which we may assume are nonzero only for the eigenvectors (indexed $\vec{v}_1$ in the above sum)), we get a solution of $Lx = y$. The terms coming from "upgrades" are part of the homogeneous solution to $Lx = 0$, and the remaining terms are some linear combination of $t^k e^{qt}, \dots, t^{k + m} e^{qt}$, as claimed.