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.