An important class of linear and time-invariant systems are modeled by linear constant-coefficient differential equations (LCCDEs), i.e. $$\frac{d^n y}{dt^n}+c_{n-1}\frac{d^{n-1} y}{dt^{n-1}}+...+c_2\frac{d^2 y}{dt^2}+c_1\frac{d y}{dt}=x(t),$$ where $y(t)$ is the output and $x(t)$ is the input.
When the input is a Dirac delta, the output is called the impulse response. If some of the coefficients $c_n$'s are complex, then it's possible for the impulse response to be complex. But if coefficients are all real, is the impulse response of such a system always purely real? If not, please provide a counterexample.
For this system to represent an LTI, you have to specify that all initial conditions must be zero. If $Y(s)$ is the Laplace transform of $y(t)$ and $x(t) = \delta(t-a)$ then we have
$$ Y(s) = \frac{e^{-as}}{s^n + c_{n-1}s^{n-1} + \cdots c_1s}. $$
Now we can factor the denominator into terms that look like $(s-r)$ for real roots or $(s-\alpha)^2 + \beta^2$ if $\alpha + i \beta$ is a complex root.
The inverse Laplace transforms of functions of the type
$$\frac{1}{s-r}$$ and
$$ \frac{1}{(s-\alpha)^2+\beta^2}$$
are well known to be real functions (see any table of inverse transforms). So we can write $Y(s)$ as a product of functions with real inverse Laplace transforms, and response $y(t)$ is a convolution of these functions, which is real-valued.