How do I find a system's impulse response from its state-space repersentation using the state transition matrix?
You can approach the problem using the state transition matrix by solving the standard non-homogeneous ODE in the first equation. The solution to $\dot{x}(t)=A x(t) + B u(t)$ is

$$x(t)=x_0 e^{At}+\int_{0}^te^{A(t-t&#39;)}Bu(t&#39;)dt&#39;$$

where $x_0=x(0)$. The quantity $e^{At}$ is called the state transition matrix (also the solution to the homogeneous ODE), which I'll refer to as $\Xi(t)$ (I don't recall the standard notation for this). Taking $x_0=0$, the equation for $y(t)$ becomes

$$y(t)=C\int_0^t\Xi(t-t&#39;)Bu(t&#39;)dt&#39;+Du(t)$$

The above equation gives you the output as the input convolved with the system impulse response and indeed, you can take the Laplace transform of the above equation to verify. Noting that the Laplace transform of $\Xi(t)=e^{At}$ is $(sI-A)^{-1}$ and convolutions in the time domain become products in the s-domain, we get 

$$Y=C(sI-A)^{-1}BU+DU$$

which gives you the same transfer function as in your question.



Regarding your comment on the fully Laplace transform approach being long, I wouldn't necessarily say it is so. However, state transition matrix approach might be simpler to implement, because several operations involving it can be computed with simple matrix multiplications and nothing more.