The way I like to do it is based on the following observation: let
$$ \bar{A} :=
\begin{bmatrix}
A & B \\
0 & 0
\end{bmatrix}, $$

where $0$ is the zero matrix (dimensions s.t. $\bar{A}$ is square). Then,
$$ \mathrm{e}^{\bar{A}t} = \begin{bmatrix} \mathrm{e}^{At} & \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B \\ 0 & I \end{bmatrix}. $$

Hence, for the integral you can just build this block matrix with $B=I$, compute the matrix exponential of it, then extract the top right block. In a more "closed" form:
$$ \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B = \begin{bmatrix}I & 0\end{bmatrix}\mathrm{e}^{\bar{A}t}\begin{bmatrix}0 \\ I\end{bmatrix}. $$

The advantage of this method with respect to using matrix inversion and/or Jordan form is that this method is numerically stable even when $A$ is singular (or close to singular). The disadvantage, of course, is that it takes a 4x bigger matrix as an input.

### Why it works

It follows from this observation: if you have the non-homogeneous ODE
$$ \dot{X}(t) = AX(t) + B, $$
its solution is
$$ X(t) = \mathrm{e}^{At}X(0) + \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B. $$

Define the auxiliary variable $U(t)$, which is constant (i.e., $U(t) = U(0)$ for all positive $t$). Then $\dot{U}(t) = 0$ and we have the system of ODEs
\begin{align*}
\dot{X}(t) &= AX(t) + BU(t), \\
\dot{U}(t) &= 0,
\end{align*}

which is a homogeneous ODE on the augmented variable $\begin{bmatrix} X(t) \\ U(t) \end{bmatrix}$. Therefore, we have

$$\begin{bmatrix} \dot{X}(t) \\ \dot{U}(t) \end{bmatrix} = \begin{bmatrix}
A & B \\
0 & 0
\end{bmatrix}\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \bar{A}\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix}$$
whose solution is
$$\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \mathrm{e}^{\bar{A}t}\begin{bmatrix} {X}(0) \\ {U}(0) \end{bmatrix},$$

but also,

$$\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \begin{bmatrix} \mathrm{e}^{At}X(0) + \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau BU(0) \\ U(0) \end{bmatrix} = \begin{bmatrix} \mathrm{e}^{At} & \int_0^t\mathrm{e}^{A\tau}\mathrm{d}\tau B \\ 0 & I \end{bmatrix}\begin{bmatrix} X(0) \\ U(0) \end{bmatrix}.$$