$e^{mx}$ in solving second order differential equations

In a book I am reading on differential equations, the author writes about the solution to a homogenous, linear, second order differential equation with constant coefficients. The author says something like, "Let us suppose that the solution is of the form $y=e^{mx}$ " . After this the author introduces the characteristic equation of a differential equation of the form mentioned above, and proceeds to describe how to solve it w/ undetermined coefficients, variation of parameters, etc.

How did mathematicians first come up with this "assumption" that the solution was of the form $y=e^{mx}$? And, are there any other forms of solutions for these types of equations?

I like to think of it as an extension of the solution of the differential equation

$$ \frac{\mathrm{d}y}{\mathrm{d}x}=ky $$

which, as can be easily checked using separation of variables, is

$$ y=y_0\mathrm{e}^{kx}. $$

Now if we have a second order ODE $$ \frac{\mathrm{d}^2y}{\mathrm{d}x^2}+a\frac{\mathrm{d}y}{\mathrm{d}x}+by=0 $$ where $a,b\in\mathbb{R}$, it is possible to re-write this as a pair of first-order ODEs. If we set $z=\frac{\mathrm{d}y}{\mathrm{d}x}$ (one of our equations) then we have the system of coupled ODEs

$$\cases{\dfrac{\mathrm{d}y}{\mathrm{d}x}=z \\\dfrac{\mathrm{d}z}{\mathrm{d}x}=-by-az} $$

and this can be solved by solving the equation $\dfrac{\mathrm{d}\boldsymbol{y}}{\mathrm{d}x}=A\boldsymbol{y}$, where $\boldsymbol{y}=(y(x),z(x))^\top$ and $$A=\left(\begin{matrix}0 & 1\\ -b & -a \end{matrix}\right)$$

and since this is a linear ODE (in more than one variable), it can be solved in an analogous way (using an integrating factor method)

\begin{aligned}\dfrac{{\rm d}\boldsymbol{y}}{{\rm d}x}-A\boldsymbol{y} & =0,\\ {\rm e}^{-\int A{\rm d}x}\dfrac{{\rm d}\boldsymbol{y}}{{\rm d}x}-A\boldsymbol{y}{\rm e}^{-\int A{\rm d}x} & =0,\\ \dfrac{{\rm d}}{{\rm d}x}\left(\boldsymbol{y}{\rm e}^{-\int A{\rm d}x}\right) & =0,\\ \dfrac{{\rm d}}{{\rm d}x}\left(\boldsymbol{y}{\rm e}^{-Ax}\right) & =0,\\ \boldsymbol{y}{\rm e}^{-Ax} & =\boldsymbol{C},\\ \therefore\boldsymbol{y}(x) & =\boldsymbol{C}{\rm e}^{Ax}. \end{aligned}

I'm sure with some more effort, it will be possible to recover linear combinations of exponentials in the individual components of the solution. Probably seems a bit rambly, but this is one way I like to convince myself where the exponential solutions come from - it is the solution to an aforementioned first-order differential equation, so by linearity, more exponentials will appear in the solution. (Don't know if that explains things well enough.)

Probably the idea comes from the first order equation, where $$ y'+ky=0 $$ could be solved by separation of variables $$ \frac{y'}{y}=-k\implies \ln y=-kx+c\implies y=Ce^{-kx}. $$