history of the double root solution of $ay''+by'+cy=0$ Motivation: It is a well-known fact that $ay''+by'+cy=0$ has solutions which are found from substituting the ansatz $y=e^{\lambda t}$ into the DEqn. It turns out that we replace the calculus problem $ay''+by'+cy=0$ with the algebra problem of solving the characteristic equation $a\lambda^2+b\lambda+c=0$. When the solution is a conjugate pair of complex numbers or distinct pair of real numbers the solutions arise from $e^{\lambda t}$. On the other hand, when the solution is real and repeated then the ansatz solution $y=e^{\lambda t}$ only covers half of the general solution.
Suppose that $a\lambda^2+b\lambda+c=0$ has double root solution $\lambda = r$ then we form the general solution of $ay''+by'+cy=0$ as
$$ y(t) = c_1e^{rt}+c_2te^{rt}. $$
The inclusion of the $t$ in the solution is surprising to many students. I think many have asked "where'd the $t$ come from?". Of course, we could just as well ask "where the $e^{\lambda t}$ come from?". I know of several ways to derive the $t$. In particular:


*

*$y''=0$ integrates twice to $y=c_1+tc_2$ and $e^{0t}=1$ so this is an example of the double root. A simple change of coordinates allows this derivation to be extended to an arbitrary double-root.

*reduction of order to a system of ODEs in normal form. We'll obtain a $2 \times 2$ matrix which is not diagonalizable. However, the matrix exponential gives a solution and the generalized e-vector piece generates the $t$ in the second solution.

*you can use the second linearly independent solution formula from the theory of ODEs. This formula is found by making a reduction of order based on the fact $y=e^{rt}$ is a solution. After a bit the problem reduces to a linear ODE which integrates to give a lovely formula with nested integrals. This formula also will derive the $t$ in the double root solution.

*Laplace transforms. We can transform the given ODE in $t$ to obtain an algebra equation with $(s-r)^2Y$ which gives $\frac{F(s)}{(s-r)^2}$ and upon inverse transform the appearance of the $(s-r)^2$ in the denominator gives us the $te^{rt}$ solution

*Inverse operators. By writing the given ODE as $(D-r)^2[y]=0$ we can integrate in a certain way and again derive the $te^{rt}$ solution.

*Series solution techniques.

*added 10/6: start with the distinct root solution $y=c_1e^{\lambda_1 t}+c_2e^{\lambda_2t}$ and consider the limit $\lambda_1 \rightarrow \lambda_2$ to derive the second solution.
These are the methods which seem fairly obvious in view of the introductory course (up to notation, several of these are the same method).  My question is this:

Question: What is the history of the solution $y=te^{rt}$? Who studied the problem $ay''+by'+cy=0$ and found this solution? 

I'm also interested in the particular sub-histories of the other methods I mention above. 
Thanks in advance for any insights!
 A: From the extraordinary book of E. L. Ince Ordinary Differential Equations, Chapter V (1956):

5·21. Fundamental Sets of Solutions.—Any linearly independent set of $n$ solutions $u_1,\,u_2\,...,\,u_n$ of the equation
  $$
L(u) = 0
$$
  is said to form a fundamental set or fundamental system*. Conversely, the condition that any given set of $n$ solutions should be a fundamental set is that the Wronskian of the $n$ solutions is not zero. The general solution of the equation will be†
  $$ 
u = C_1 u_1 + C_2 u_2 + ... + C_n u_n,
$$
  which cannot vanish identically unless the constants $C_1,\,C_2\,\,...C_n$ be all zero.
  $${}$$
  $^*$ The term fundamental system is due to Fuchs, J. für Math. 66 (1866), p. 126
† Lagrange, Misc. Taur., 3 (1762-65), p. 181

And goes on constructing a fundamental set and proving -very formally I might add- that satisfies the differential equation. Then follows

When the fundamental system is taken in the form
  $$
u_1 = v_1,\, u_2 = v_1 \textstyle{\int} v_2 dx,\, ...,\, u_n = v_1  \textstyle{\int} v_2 \textstyle{\int} ... \textstyle{\int}v_n (dx)^{n-1},
$$
  the equation becomes‡
  $$
\frac{d}{dx}\cdot\frac{d}{v_n dx} \cdot \frac{d}{v_{n-1} dx} \cdots \frac{d}{v_2 dx} \cdot \frac{u}{v_1} = 0.
$$
  $${}$$
  ‡ Frobenius, J. für Math. 76 (1873), p.264; 77 (1874), p. 256.

Ending the section proving that the differential operators are not in general permutable. Then

5·22. Depression of the Order of an Equation.—If $r$ indepentent solutions of the equation of order $n$,
  $$
L(u) = 0,
$$
  are known, then the order of the equation may be reduced to $n-r$. For let 
  $$
u_1,\, u_2\, ...,\, u_r
$$
  be the known solutions, and let
  $$
v_1 = u,\, v_2 = \frac{d}{dx}\left(\frac{v_1}{u_2}\right),\, v_3 = \frac{d}{dx}\left\{\frac{1}{v_2}\cdot\frac{d}{dx}\left(\frac{u_3}{v_1}\right)\right\},
$$
  and so on as before. Then since the equation is known to be ultimately of the form
  $$
\frac{d}{dx} \cdot \frac{d}{v_n dx} \cdot ... \cdot \frac{d}{v_{r+1} dx} \cdot \frac{d}{v_r dx} \cdot ... \cdot \frac{d}{v_2 dx} \cdot \frac{u}{v_1} = 0,
$$
  it might be written as
  $$
P(v) = 0,
$$
  where
  $$
v = \frac{d}{v_r dx} \cdot \frac{d}{v_{r-1} dx} \cdot ... \cdot \frac{d}{v_2 dx} \cdot \frac{u}{v_1} \tag{A},
$$
  and $P$ is a linear operator of order $n-r$.
If any solution of $P(v) = 0$ is obtainable, the corresponding value of $u$ may be obtained directly from (A) by $r$ quadratures.

Then

5·23. Solution of the Non-homogeneous Equation.—Consider now the general equation
  $$
L(x) = r(x), \tag{A}
$$
  it being supposed that a fundamental set of solutions $u_1(x),\,u_2(x),\, ...,\, u_n(x)$ of the reduced equation are known.
Then the general solution of the reduced equation is
  $$ 
u = C_1 u_1 + C_2 u_2 + ... + C_n u_n,
$$
  in which $C_1,\, C_2,\, ...,\, C_n$ are arbitrary constants. Now, just as in the case of linear equations of the first order (§2·13), so here also the method of variation of parameters§ can be applied to determine the general solution of the equation under consideration.
  $${}$$
  § Lagrange, Œuvres, 4, pp. 9, 159

In Chapter IV:

6·1. The linear Operator with Constant Coefficients.—The homogeneous linear differential equation with constant coefficients
  $$
A_0 \frac{d^n y}{dx^n} + A_1\frac{d^{n-1} y}{dx^{n-1}} + ... + A_{n-1} \frac{d y}{dx} + A_n y = 0 \tag{A}
$$
  was the first equation of a general type to be completely solved¶.
  $${}$$
¶ It appears that the solution was known to Euler and to Daniel Bernoulli about the year 1739. The first published account was given by Euler, Misc. Berol. 7 (1743), p. 193; see also Inst. Calc. Int. 2, p.375

Next

6·14. The Case of Repeated Factors regarded as a Limiting Case.—A very powerfull method of attacking the case in which the operator
  $$
A_0 D^n + A_1 D^{n-1} + ... + A_{n-1} D + A_n
$$
  has a repeated factor is due to d'Alembert**.
  $${}$$
  $^{**}$ Hist. Acad. Berlin 1748, p. 283.

How about that for a little history!
(The reproduction of the sections was intended to be ilustrative and for educational purposes only. No copyright violation was intended. Please, only edit if you find typos or grammar. The obscure notation is part of the text and reflects the historical trend)
A: Here are some relevant pieces of history:
Euler wrote to Johann Bernoulli in 1739, describing how to solve second or higher order linear differential equations with constant coefficients, such as
$$ y+a{dy\over dx}+b{d^2y\over dx^2}+c{d^3y\over dx^3}+\dots=0 $$ 
He did this by introducing the characteristic equation, and factoring it into linear and quadratic terms. In slightly modernized notation, this is given by
$$ 1-ap+bp^2-cp^3+\dots=\prod_i(1-\gamma_ip)\cdot \prod_j (1-\alpha_j p+\beta_j p^2)$$ 
The linear factors produce solutions on the form $C{\rm e}^{-x/\gamma_i}$, and the quadratic factors produce solutions involving trigonometric functions. The case with multiple roots is not mentioned in this letter.
This letter is available online, with the relevant part beginning on p. 37. It is a highly interesting and very important letter, as it brings forth the now familiar connection between exponential and trigonometric functions, and it also establishes sine and cosine as functions rather than line segments in geometric figures.
${}$
An explicit reference is given by Euler's Integral Calculus, published in three volumes (1768-1770). The third edition is available online, and in the Second Part, Chapter 4, p. 79, Problem 102, he considers explicitly the case with multiple roots. He says that the general solution of the equation
$$\partial \partial y + A\partial y \partial x + By \partial x^2 =0$$
is given by 
$${\rm e}^{-{1\over 2}Ax}\bigl(\alpha{\rm e}^{nx}+\beta{\rm e}^{-nx}\bigr)$$
where $n=\sqrt{A^2/4-B}$, and that if $n=0$ then this produces the result
$${\rm e}^{-{1\over 2}Ax}\bigl(\alpha+\beta x\bigr)$$
${}$
A date in the 1760's must be much too late for the discovery, however. In fact, according to Victor Katz, p. 554, Euler was familiar with the exponential function already by 1730, and he used it in solving various differential equations in the 1730's. Quite possibly equations such as $y''+2y'+y=0$ were solved then.
A: Years ago when I was looking at a related problem, I took the limit as the 2 roots approached each other. That led to the pair of solutions.
A: In Arnold's Ordinary Differential Equations there is a passage about history (Section 26.4 in the latest Russian edition). It says something like this. If we have two roots $\lambda_1$ and $\lambda_2$ then geometrically $c_1e^{\lambda_1t}+c_2e^{\lambda_2t}$ is a plane. If $\lambda_1\to \lambda_2$ then this plane turns into a straight line. Question: is it possible to find a limiting position of the original plane? Consider another basis for the plane: $e^{\lambda_1t}$ and $e^{\lambda_1t}-e^{\lambda_2t}\approx (\lambda_2-\lambda_1)te^{\lambda_1t}$, hence when $\lambda_1\to\lambda_2$ then $(e^{\lambda_1t}-e^{\lambda_2t})/(\lambda_1-\lambda_2)\to te^{\lambda_1 t}$. That's, according to Arnold, how Euler and Lagrange guessed the solution $te^{\lambda t}$.
A: I think this is straightforward enough that any of the greats of that period could have figured it out. Consider the differential equation in the form
$$(D - r)^2y = 0$$
If you let $z = (D - r)y$, this becomes
$$(D - r)z = 0$$
So $z(x) = Ce^{rx}$, and the original differential equation becomes
$$(D - r)y = Ce^{rx}$$
This is a first-order linear equation; multiply through by $e^{-rx}$ and it becomes
$$e^{-rx}y' - re^{-rx}y = C$$
The left hand side is $(e^{-rx}y)'$ and can now be readily solved.
I look at this and I imagine someone like Euler or Bernoulli would have very quickly figured something like this out. Such elementary calculus arguments were really their specialty. So I doubt there is one person who first noticed it. I would not be surprised if Isaac Newton himself had the occasion to figure it out, since such differential equations come up so often in classical mechanics. Of course some of the other ways of solving it that have been mentioned here might have occurred to them before what I wrote here.
A: I personally like method 3.  Many books just say "try this" for the times t but I don't think it is intuitive that this would work--for instance why is it never in the non repeated root version?  But you really don't need to guess as knowing a single solution, the other can be derived.  My copy of Thomas Finney (4th edition, 1980s) does this well but many books (Tenenbaum, Granville, etc.) just use the try it, it works explication for the repeated root with times t.  
Conversely, I do think guessing with exponentials (or sin/cos) is reasonable as a start to the whole problem as the reader will be familiar with how they don't disappear (like polynomials) or reduce in order, from the many calc problems with them.  (Even integration by parts sees this often, the disappearing versus not disappearing).  Therefore a solution based on these is likely to satisfy the equation (with appropriate lambda to make the terms cancel). From there it is short work to the characteristic equation itself.
