Given some linear constant coefficient DE with:
$$ (a_0D^n+a_1D^{n-1}+...+a_{n-1}D+a_n)y=f(t) \tag{1} $$
Without loss of generality, we may assume that $a_0>0$, then: $$ \textit{(1) is stable} \Rightarrow a_0,...,a_n >0 $$
I've found the above statement on the MIT online course on ODE's. Since there is no analytic way of finding the roots of an n-order polynomial, I've tried proving it by rewriting the characteristic equation in Horner's scheme for some root $Re{\lambda_i}<0$ (since it's stable, all roots must have negative real part) - but without success. I also lack any intuition for why this statement should be true. Any hints on how to prove it?
The characteristic equation of the homogenous DE is
$$P(x) = a_0 x^n + a_{1}x^{n-1} + \ldots + a_n$$
As you say; the DE is stable if all the roots of $P(x)$ have negative real part (this is the Routh–Hurwitz stability criterion). One way to see this is as follows: if a root $\lambda$ has positive real part then the general solution of the DE will contain a term $Ae^{\lambda t}$ which grows without bounds as $t\to\infty$. Now pick initial conditions such that $A = 0$ then (assuming for simplicity here that $\lambda$ is the only root with positive real part) the solution will be bounded for all $t$. However if we tweek the initial conditions just slightly such that $A \not = 0$ then no matter how small $A$ is the solution will diverge as $t\to\infty$. This is incompatible with stability.
Now to the proof. Assume all roots of $P(x)$ have $\text{Re}(\lambda) < 0$ and $a_0 > 0$. We start by factorizing $P(x)$. If a root $\lambda$ is real then $P(x)$ is divisible by $(x-\lambda)$. If $\lambda$ is imaginary then since $P(x)$ is a real polynomial $(a_i\in \mathbb{R})$ we have that the complex conjugate $\lambda^*$ is also a root so $P(x)$ is divisible by
$$(x-\lambda)(x-\lambda^*) =x^2 - 2\text{Re}(\lambda)x + |\lambda|^2$$
It follows that we can factorize $P(x)$ as follows:
$$P(x) = a_0\prod_{i\,:\,\lambda_i \in\mathbb{R}} (x-\lambda_i)\cdot \prod_{i\,:\,\text{Im}(\lambda_i)>0}(x^2 - 2\text{Re}(\lambda_i)x + |\lambda_i|^2)$$
Since we know that $\text{Re}(\lambda_i) < 0$ we have that $x^2 - 2\text{Re}(\lambda)x + |\lambda|^2 = x^2 + 2|\text{Re}(\lambda)|x + |\lambda|^2$ so we see that all the terms in all the factors in the expression above are positive so if we now imagine multiplying out the product above we would get
$$P(x) = a_0 x^n + (\text{positive})x^{n-1} + \ldots + (\text{positive})$$
and it follows that $a_0,a_1,\ldots,a_n > 0$.