Denote the convolution $f*g(u)=\int_0^u f(\tau)g(u-\tau)d\tau$. For matrix functions $A$ and $B$, we have
$$A*B(u)=\int_0^u \sum_j A_{ij}(\tau)B_{jk}(u-\tau)d\tau=\sum_j A_{ij}*B_{jk}(u). \tag{$\circ$}$$
Also, we have the rule
$$\mathcal{L}\{f*g\}(s)=\mathcal{L}\{f\}(s)\cdot\mathcal{L}\{g\}(s)=F(s)G(s). \tag{a}$$
Inverting this, we deduce
$$\mathcal{L}^{-1}\{F(s)G(s)\}(x)=f*g(u). \tag{b}$$
Use summation notation and linearity of the operator (and remember $(\circ)$) to prove that $\rm(a)$ and $\rm(b)$ also apply to matrix functions as well. This will help in your problem if you set
$$F(s)=(sI-A)^{-1}x, \\ G(s)=y^T(sI-A)^{-1}.$$
Check with summation notation that, again by linearity,
$$\mathcal{L}\{A(u)C\}=\mathcal{L}\{A(u)\}C, \\ \mathcal{L}\{CA(u)\}=C\mathcal{L}\{A(u)\} \tag{c}$$
for matrix functions $A(u)$ and constant matrices $C$. Notice how $(c)$ can be reworked for the inverse transform just as well. This will help with the constant matrices $x$ and $y^T$.
(Throughout this answer, I have taken matrices to be of arbitrary dimensions - outside of the multiplications being well-defined in each instance.)
This doesn't actually do the work for you, but it tells you what you need to do. The final thing you'll need to do is a convolution involving (but not quite "of," due to $x$ and $y$ in the mix) matrix exponentials. Remember that if $S,T$ are commuting matrices, $e^{S+T}=e^Se^T$ is valid.
To actually evaluate the inverse Laplace transform integral, we must perform an integration in the complex plane. Many times, this means using the theorems of complex analysis, notably Cauchy's residue theorem. Because many students are introduced to the Laplace transform before they take complex analysis - many times, such students won't ever take a course in complex analysis - most first time lessons focus on using tables, to your chagrin.
If you are comfortable with integration in the complex plane, however, evaluating the ILT provides a terrific exercise for your contour integration. For example, consider a LT $\hat{f}(s)$ having poles in the complex plane. Rather then try to actually evaluate the ILT integral directly, we consider the integral over the closed contour $C$:
$$\frac{1}{i 2 \pi} \oint_C ds \, \hat{f}(s) \, e^{s t}$$
What is $C$? It turns out that our choice of $C$ will depend on the sign of $t$. When $t>0$, we will need to add a circular contour to the left of the vertical line of integration of the ILT. When $t<0$, we close to the right instead. The reason for this is that we need the integral over the circular contour to vanish as the radius $T$ of the circle approaches $\infty$. To see this, parametrize the integral there as $s = T e^{i \phi}$; then the integral over a leftward-circular contour looks like
$$i T\int_{\pi/2}^{3 \pi/2} d\phi \, e^{i \phi} \, \hat{f}(T e^{i \phi}) e^{t T \cos{\phi}} e^{i t T \sin{\phi}}$$
Note that $\cos{\phi} < 0$ within the integration region; thus the integral only vanishes when a) $t > 0$, and b) $\lim_{T \to \infty} T \hat{f}(T) = 0$.
This brings us to my next point: causality. The ILT $f(t)$ must be zero when $t<0$. This means that $\gamma$ must be chosen so as to be to the right of the pole with the largest real part. That way, when the contour closes to the right, the integral over the closed contour will be zero (as it encloses no poles) and, hence, the ILT is zero by design. Note that it doesn't matter what $\gamma$ is, so long as there are no poles to the right of the line $\Re{z}=\gamma$.
In the simple cases we find from, for example, solving ordinary differential equations with constant coefficients, we will have $\hat{f}(s)$ being a rational function of $s$. In such cases, we will simply have
$$\frac{1}{i 2 \pi} \int_{\gamma-i \infty}^{\gamma+I \infty} ds \, \hat{f}(s) \, e^{s t} = \sum_k \text{Res}_{s=s_k} [\hat{f}(s) \, e^{s t}]$$
As long as you know how to evaluate these residues, then evaluating such ILTs is straightforward. Let's look at a concrete example: solve the differential equation
$$x'' + 4 x= 1$$
with $x(0) = x'(0) = 0$. We apply LTs to both sides and we get
$$\hat{x}(s) = \frac{1}{s (s^2+4)}$$
To obtain $x(t)$, we sum over the residues of $\hat{x}(s) e^{s t}$ at the poles at $s=0$, $s=2 i$, and $s=-2 i$. At $s=0$, the residue is
$$\lim_{s \to 0} s \frac{e^{s t}}{s (s^2+4)} = \frac14$$
The other residues are found similarly:
$$\lim_{s \to 2 i} (s-2 i) \frac{e^{s t}}{s (s^2+4)} = \frac{e^{i 2 t}}{-8}$$
$$\lim_{s \to -2 i} (s+2 i) \frac{e^{s t}}{s (s^2+4)} = \frac{e^{-i 2 t}}{-8}$$
The ILT is the sum of the residues, or
$$x(t) = \frac{1-\cos{2 t}}{4} = \frac12 \sin^2{t}$$
There are times, though, when you need other techniques of contour integration. For example, sometimes you may encounter a LT having a branch point in $s$. One example I have encountered here is $\hat{f}(s) = e^{-\sqrt{s}}$. Here is my evaluation of this ILT. Note how I treated the contour. Despite the fact that there are no poles and hence no residues, the multivaluedness of the LT provides the ILT I sought. I have many more worked examples here.
I hope this helps. Please feel free to ask questions about this topic anytime.
Best Answer
Intuitively, the derivative of the Dirac delta function $\delta'$ has Laplace transform $s$. The derivative of the Dirac delta is a generalized function that pulls out the derivative of the function with a change of sign: for any interval $[a,b]$ where $a < 0 < b$,
$$\int_a^b \delta'(t)f(t) \ dt = \left[ \delta(t) f(t) \right]_a^b - \int_a^b \delta(t) f'(t) \ dt = -f'(0)$$
Applying that procedure inductively,
$$L^{-1}\{s^n\}(t) = \delta^{(n)}(t)$$