Where does $xe^x$ solution come from when the characteristic polynomial is square? When solving the differential equation $y'' + ay' + by = 0$ (with constant, real coefficients $a$ and $b$, although they could be complex if you like), you do it by setting up the characteristic equation $r^2 + ar + b = 0$, finding its solutions $r_1, r_2$, and then the general solution to this equation is $Ce^{r_1x} + De^{r_2x}$. This works both when the solutions are real and when they are complex.
However, when we have a double root $r_1 = r_2$, we get a different general solution, namely $Ce^{r_1x} + Dxe^{r_1x}$. I have no trouble seing that this is indeed a solution, and intuitive reasoning on degrees of freedom dictates that we must have a linear combination of two terms in our general solution, while $e^{r_1x}$ and $e^{r_2x}$ are the same. So the fact that there is a second term of some other form is not surprising.
I have, however, yet to see a "natural" explanation of this $xe^{r_1x}$ term. If one were developing the theory from scratch, how would one find this solution (other than blind luck)? If I wanted to teach ODE's to a class of students "the right way", i.e. with good explanations and motivations for everything (as opposed to just pulling out ready-made solutions like what was done to me when I was learning this exact thing), how would I motivate even considering a term like $xe^{r_1x}$ (other than "Well, exponentials aren't quite cutting it, but this is kindof like an exponential, right? Let's try it.")? And is there a way of solving the general differential equation that does not involve splitting into cases depending on whether the characteristic polynomial is a square?
 A: Here is an algebraic way to do this. Suppose that the caracteristic equation of $y'' + ay' + y =0$ has a double root $r$. This means that $X^2 + aX + b = (X-r)^2$, hence $y'' + ay' + y = A^2(y)$, where $A$ is the endomorphism :
$$ A(y) = (D-r)(y) = y' - ry,$$
where $D$ is the usual derivation.
The kernel of $A$ is well known. Thus $A(e^{rx})=0$.
Note that $A$ satisfies the Leibnitz rule in the following sense :
$$ A(fg) = f' \times g + f \times A(g).$$
Hence in order to solve $A^2(y)=0$, write $y(x) = z(x)e^{rx}$. The Leibnitz rule and the fact that $A(e^{rx})=0$ imply :
$$ A(y) = z'(x) e^{rx}.$$
$$ A^2(y) = z''(x) e^{rx}.$$
So $A^2(y)=0$ iff $z''(x)=0$.
Some highlights : there is the notion of a differential module $M$ over a differential ring $R$. Here the ring $R$ is $C^\infty(\mathbb{R})$ (or you can also take $\mathbb{R}[X]$) with usual derivation and the module $M$ is $C^\infty(\mathbb{R})$ where the derivation is $A=D-r$. Here I have just found a basis in order to make $M$ isomorphic to $R^n$.
A: I don't know if you think it's "natural", but the $x$ comes from the method of reduction of order.  If $y_1$ is a one solution of a linear, homogeneous equation and we need a second, linearly independent solution, a reasonable guess is $y_2 = v(x)y_1$.  Sort of as you explain it, $ce^{rx}$ is a solution for all constants $c$, but they're all linearly dependent.  So we keep the "solution-ness" of the $e^{rx}$ but remove the "constant-ness" of the $c$ by replacing it with a function $v(x)$.  Then by straightforward calculations, we find, in your case, that $v=x$ works.
A: The simple direct approach gives the solution $xe^x$ without much hassle. Let the equation be $$y''-2y'+y=0$$ and let $z=y'-y$ so that the equation can be written as $$z'-z=0$$ The above equation on multiplying with $e^{-x} $ gives $$(ze^{-x}) '=0$$ or $$ze^{-x} =c_1$$ so that $$y' - y=z=c_1e^x$$ Again multiplying by $e^{-x} $ gives us $$(ye^{-x}) '=c_1$$ so that $$ye^{-x} =c_1x+c_2$$ or $$y=c_1xe^x+c_2e^x$$
A: I have also wondered about this question for some time, but the most satisfactory answer that I got was that the following:
if $r_0$ is a double root of the characteristic equation of a differential equation of constant coefficient, i.e 
$$L(e^{rx}) = p(r) e^{rx},$$
where $p$ is the characteristic polynomial of the ODE correspond to $L$, and $p(r_0) = 0$. Then $p$ has the form 
$$p(r) = A \cdot (r-r_0)^2.$$
Now observe that 
$$p'(r) = A \cdot 2 \cdot (r-r_0),$$
hence $r_0$ is also a root of $p'(r)$.
Therefore, if we differentiate $L$ wrt $r$, we get
$$\frac{dL(e^{rx})}{dr} = \frac{dL(e^{rx})}{ d(e^{rx}) } * \frac{d((e^{rx}))}{dr } = p(r) \cdot xe^{rx}  = L(x \cdot e^{rx})$$
and also
$$\frac{dL(e^{rx})}{dr}= \frac{d(p(r)e^{rx})}{dr} = [p'(r) + xp(r)]\cdot e^{rx},$$
hence
$$L(x \cdot e^{rx}) = [p'(r) + xp(r)]\cdot e^{rx},$$
and plugging $r = r_0$, we see that 
$$L(x e^{r_0 x}) = 0,$$
since both $p'(r_0)$ and $p(r_0)$ equal to zero.
Therefore, $x e^{r_0 x}$ is also a solution of the ODE given by $L$.Moreover, since $e^{rx}$ and $xe^{rx}$ are linearly independent, we have two independent solution from a double root.

tl:dr
When the the characteristic equation has double root $r_0$, the very fact that $r_0$ is both a zero of $p(r)$ and $p'(r)$ makes $xe^{r_0 x}$ another independent solution.
A: L'Hôpital's rule. 
If you take a solution, with given initial conditions, for a differential equation where the two roots of the characteristic polynomial differ, and let one root approach the other, then the limit (using L'Hôpital) has that term like $xe^x$ in it.
Example.  The differential equation
$$
y'' - (1+a)y'+ay=0,\quad y(0)=0, y'(0)=1\tag{$$1}
$$
with $a \ne 1$ has characteristic equation $r^2-(1+a)r+a$ with zeros $1,a$.  The solution is
$$
y = \frac{-e^x+e^{ax}}{a-1}\tag{$2$}
$$
Now set $a=1$ in ($1$) ... the differential equation
$$
y'' - 2y'+y=0,\quad y(0)=0, y'(0)=1\tag{$1'$}
$$
has characteristic equation $r^2-2r+1$ with zeros $1,1$.  The solution is
$$
y = xe^x\tag{$2'$}
$$
Note that, using L'Hôpital's rule, the limit of ($2$) is ($2'$).
$$
\lim_{a \to 1} \frac{-e^x+e^{ax}}{a-1} = xe^x .
$$
A: This relies on much of the same intuition already given in user10676's excellent answer. However, I hope that this is different enough to have its own independent merits.
A homogeneous linear ODE with constant coefficients can be represented as a series of operations of the form $ (D-r) $, wherein $ D $ represents differentiation and $ -r $ represents ordinary multiplication by the negative of the root. For example, if a constant-coefficient homogeneous linear ODE has roots 1, 1, and 3 it can be represented as
$$ (D-3)(D-1)(D-1)y(x)=0 $$
In order for $ y(x) $ to be a solution to the differential equation with the double root $ r $, it must satisfy the equation:
$$
(D-r)(D-r)y(x)=0\tag{$1$}
$$
$ e^{rx} $, the solution when the operator $ (D-r) $ is applied only once (when the root is not repeated), is already a solution because $ (D-r)e^{rx}=0 $ and $ (D-r)0=0$. However, we may find another solution by finding the solution to the equation:
$$ (D-r)y(x)=e^{rx}\tag{$2$} $$
because the $ (D-r) $ operator applied another time to its output will also yield 0. In other words, if (2) applies, $ e^{rx} $ can be substituted for $ (D-r)y(x) $ like so:
$$ (D-r)(D-r)y(x)=0 $$
$$ (D-r)e^{rx}=0 $$
which is already known to be true.
In order to find what function satisfies (2), some guesswork is employed. This does make this deduction not fully straightforward, but I hope the negative effect on how "natural" this seems isn't especially large. By a simple rearrangement, we see that the following must be true:
$$ \frac{dy(x)}{dx} = e^{rx} + ry(x) $$
If we make the educated guess that this is the derivative of the product of two functions, $ a(x) $ and $ b(x) $, according to the Leibniz rule, then it can be assumed from the term $ ry(x) $ that
$$ \frac{db(x)}{dx}=rb(x) $$
$ b(x) $ is recognizable as the function $ e^{rx} $. Therefore, our function will take the form of $ a(x)e^{rx} $. Furthermore, we can now deduce from the first term of the derivative, $ e^{rx} $, that
$$ \frac{da(x)}{dx}=1 $$
which is satisfied when $ a(x)=x $. Therefore our solution is
$$ y(x)=a(x)b(x)=xe^{rx} \tag*{$\square$} $$
The $ (D-r) $ operators for roots other than the double root don't interfere with the $ xe^{rx} $ as a solution, as, for example:
$$ (D-r_1)xe^{rx} = e^{rx} + rxe^{rx} - r_1xe^{rx} $$
which would also constitute a solution to (1).
A similar approach works for roots repeated more than twice as well, as $ (D-r)x^2e^{rx} = 0.5xe^{rx} $ and so on. As one increments the number of times the root is repeated, a new term with a one-higher power of $ x $ becomes a solution as well.
