Splitting Numerical Schemes

Suppose we have a numerical scheme

for some matrix .
Often, is naturally represented as .
In Numerical Schemes as Approximations to Exp,
we want to find
Using Product of matrix exponentials:

and if and commute:

Inhomogeneous case

The above works for zero boundary conditions.
If that was not the case, we would have something like:

where originates from the boundary conditions.
The exact solution for this is:

we can approximate this by:

Then for example use (integral) trapezoidal rule to find:

and then use the splitting.