desolve function: how is the matrix inverse of the Laplace transform computed?

Hi all,

I'm using the desolve function to solve a system of linear differential
equations in the form of y'(t) = A*y(t) + b(t).
If I understand Maxima's documentation well, the desolve function uses
the Laplace transform so that the above system is seen as sY(s) - y(0) =
A*Y(s) + B(s) where Y(s) and B(s) are the Laplace transforms of y(t) and
b(s) respectively.
Thus, desolve should compute: Y(s) = (sI-A)^(-1) * (B(s) + y(0))

My question is: how is (sI-A)^(-1) computed? Does it rely on the
analytic formula based on matrix cofactors? If so, how can the roots of
the determinant's polynom be algebraically computed for any number of
equations? How well does it scale to large systems (let say, more than
100 equations) ?

Otherwise, if the computation does not rely on the computation of matrix
cofactors, what is the method used by desolve? Are there any academic
papers that explain the method?

Thank you in advance for your explanations. Please receive my apologies
if the method is already described or referenced in a place that I was
not able to find.

Florent Teichteil.