Linear homogeneous recurrence relations with repeated roots; motivation behind looking for solutions of the form $nx^n$?
Solution 1:
Linear Algebra:
To set the stage, let $V$ be the vector space of sequences $a_0, a_1, ... $ (say of elements of an algebraically closed field) indexed by non-negative integers. On this vector space define the shift operator $S : V \to V$, which sends a sequence $a_i$ to the sequence $(Sa)_i = a_{i+1}$. To say that $a$ satisfies a linear recurrence with characteristic polynomial $p$ is then precisely to say that $$p(S) a = 0.$$
Since we're working over an algebraically closed field, we can factor $$p(S) = \prod_{i=1}^n (S - \lambda_i I)^{m_i}$$
and then write our recurrence as $$\prod_{i=1}^n (S - \lambda_i I)^{m_i} a = 0.$$
In particular, any sequence satisfying $(S - \lambda_i I)^{m_i} a = 0$ (a generalized eigenvector of $S$ with eigenvector $\lambda_i$) satisfies the recurrence, and standard results in linear algebra imply that generalized eigenvectors span the entire space of solutions.
So it suffices to solve the generalized eigenvector equation. First, I'm not sure if you know this, so it's worth saying: if $m_i = 1$, then the space of solutions to $$(S - \lambda_i I) a = 0$$
is spanned by $a_n = \lambda_i^n$ and this is the clearest reason I can think of (eigenvector decomposition) to look for exponential solutions.
Before we tackle the general case, another special case: $\lambda_i = 1$. In this case $S - I$ is just the forward difference operator, and the sequences satisfying $$(S - I)^n a = 0$$
are precisely the polynomial sequences of degree less than $n$. This should be a familiar fact about the forward difference operator, and it turns out to imply the general result. Indeed, suppose $a$ satisfies $$(S - \lambda_i I)^n a = 0$$
and write $a_n = \lambda_i^n b_n$ for some sequence $b_n$.
Lemma: $((S - \lambda_i I)a)_n = \lambda_i^{n+1} ((S - I)b)_n$.
Proof. Straightforward computation.
By induction, it follows that $$(S - \lambda_i I)^n a = 0 \Leftrightarrow (S - I)^n b = 0$$
and the conclusion follows. (I haven't discussed the case $\lambda_i = 0$ but the solutions are easy to describe here.)
A Limiting Argument:
We'll work through a simple representative example since the goal here is motivation. Consider a sequence $a$ satisfying $(S - \lambda_1 I)(S - \lambda_2 I) a = 0$ where $\lambda_1 \neq \lambda_2$, and let's see what happens if we let $\lambda_1 \to \lambda_2$. To take this limit sensibly we'll fix initial conditions $a_0 = 0, a_1 = 1$. This will give $$a_n = \frac{\lambda_1^n - \lambda_2^n}{\lambda_1 - \lambda_2}$$
Now taking the limit $\lambda_2 \to \lambda_1$ and applying, for example, l'Hopital's rule, we obtain $$a_n = n \lambda_1^{n-1}.$$
Differential Equations:
The shift operator $S$ defined above has another incarnation as differentiation $D$: indeed, we just need to associate to a sequence $a_0, a_1, ...$ the exponential generating function $$A(t) = \sum_{n=0}^{\infty} a_n \frac{t^n}{n!}.$$
Solutions to $p(S) a = 0$ then become solutions to differential equations $p(D) A = 0$. But now that we can work in the formalism of differential equations we can appeal to physical intuition: in physical terms, the relevant phenomenon is resonance. More concretely, let's work with a forced harmonic oscillator $$m \frac{d^2 x}{dt^2} + kx = \sin \omega t.$$
(Think of periodically pushing a swing.) Note that $x$ satisfies $(m D^2 + k)(D^2 + \omega^2) x = 0$. This polynomial has repeated roots when $\omega^2 = \frac{m}{k}$, that is, when the frequency of the forcing function matches the "natural frequency" of the system (if physicists have better terms for these I've forgotten them). In that case one can observe experimentally in the relevant physical systems (e.g. the Tacoma Narrows bridge) that the amplitude of the oscillations increase linearly with time, and one can observe mathematically that the solution involves terms like $t \sin \omega t$ (and one simple way to see this is to use the limiting argument above). Computing the Taylor series of such solutions then gives the familiar extra factor of $n$.
Solution 2:
One way to look at this would be using generating functions.
The generating function comes out to
$$F(x) = \frac{Q(x)}{P(x)}$$
where the polynomial $P(x)$ is got from the characteristic polynomial, by replacing $x$ by $\frac{1}{x}$ and multiplying by an appropriate power of $x$.
Now using the partial fraction decomposition
$$ \frac{Q(x)}{(x-a_1)^{n_1} \dots (x-a_k)^{n_k}} = \sum_{i=1}^{k} \sum_{j=1}^{n_i} \frac{A_{ij}}{(x-a_1)^j}$$
tells us that we should get a polynomial in $n$.