history of the double root solution of $ay''+by'+cy=0$
Solution 1:
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)
Solution 2:
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.
Solution 3:
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}$.
Solution 4:
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.
Solution 5:
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.