Solving a first-order inhomogeneous matrix differential equation? Given non-commuting matrices $A$ and $B$ of order $n$, is there a closed-form solution to the differential equation $$\frac{dX}{dt} = AX + tBX$$ with $X(0) = I$?
I know that for the reals, $x = a \exp{\int f(t)}$ is the general solution to $\dot{x} = f(t)x$, but I'm also 99% certain this relies on the commutivity of the reals. 
I'm more specifically looking to numerically compute $X(T)$ given the more general differential equation $$\frac{dX}{dt} = f(t)X,\;\;X(0)=I$$ but in circumstances where $f'(t)$ may be large and a $1$st order piecewise approximation would be far more accurate than $0$th order for any given $\Delta t$. Ultimately my concern is computing $X(T)$ as quickly as possible.
Are there better techniques for accomplishing this?
 A: You may use the Dyson series to solve the equation like
$$\frac{d}{dt}X(t)=A(t)X(t)$$
Where $X:\mathbb{R}^{+}\rightarrow\mathbb{C}^{n}$ and $A:\mathbb{R}^{+}\rightarrow\mathbb{M}_{n\times{n}}(\mathbb{C})$. First let 
$$X(t)=U(t)X(0)=U(t)I$$
Where $U:\mathbb{R}^{+}\rightarrow\mathbb{M}_{n\times{n}}(\mathbb{C})$, such that $U(0)=I$ (to match the initial conditions), thus
$$\frac{d}{dt}U(t)=A(t)U(t)$$
This may be converted to an integral equation
$$U(t)=I+\int_{0}^{t}dt_{1}A(t_{1})U(t_{1})$$
The von-Neumann expansion than may be written
$$U(t)=I+\int_{0}^{t}dt_{1}A(t_{1})+\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}A(t_{1})A(t_{2})+\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}A(t_{1})A(t_{2})A(t_{3})+..$$
$$=I+\sum_{k=1}^{\infty}\int_{0}^{t}dt_{1}...\int_{0}^{t_{k-1}}dt_{k}\mathcal{T}[A(t_{1})...A(t_{k})]$$
This series is known as the dyson series. There $\mathcal{T}[.]$ is the time ordering operator, i.e. $\mathcal{T}[A(t_{3})A(t_{1})A(t_{2})]=A(t_{1})A(t_{2})A(t_{3})$, it is nessesary as the matricies taken at different times do not commute, e.g. in your example $[A+tB, A+t'B]=(t-t')[B, A]$. 
The convergence condition for $t\in[0, \tau]$ may be shown to be
$$\int_{0}^{\tau}||A(t)||_{2}dt<\pi$$
One may additionally show that this series is equaivalent to the so called time ordered exponential
$$U(t)=\mathcal{T}\Big[\exp\Big(\int_{0}^{t}A(s)ds\Big)\Big]=\exp\Big(\sum_{k=1}^{\infty}\Omega(t)\Big)$$
With
$$\Omega_{1}(t)=\int_{0}^{t}dt_{1}A(t_{1})$$
$$\Omega_{2}(t)=\frac{1}{2}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}[A(t_{1}), A(t_{2})]$$
$$\Omega_{3}(t)=\frac{1}{6}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}([A(t_{1}), [A(t_{2}), A(t_{3})]]-[A(t_{3}), [A(t_{2}), A(t_{1})]])$$
$$...$$
This series is known as Magnus series, has the same convergence conditions and may be shown to be equivalent to the Dyson one.
