Solving $X'(x)=X(x)M(x)$ I am aware of how to solve a system $x'=Ax$ for an $n$-vector $x$ and $n\times n$ matrix M and the necessary generalization for matrices. If $X \in \mathbb{R}^{2\times2}$ is a matrix, forming $X'(x) = M(x)X(x)$.
Question: How would I go about solving a system of the form $X'(x) = X(x) M(x)$?
Edit: An equation $X(x) = C(e^{\int M(x) dx})^{-1}$ will not suffice as $X'(x)\neq -X'(x) M(x)$ due to lack of commutativity.
 A: There is an explicit series formula for the solution generalizing the exponential series for the case where $M$ is constant. Once you have the right formula it is quite easy to check it works, but how to get the right formula in the first place? Let's start by writing
$$ X(x) = \sum_{n=0}^\infty X_n(x). $$
where $X_n(x)$ are to be determined.
You didn't specify an initial condition, but a natural one is $X(0)=I$, the identity matrix (the solutions for other initial conditions can be easily obtained from the solution for this, one by the way). A simple way to make the initial condition hold from our Ansatz is to take
\begin{align*}
X_0 \equiv I && \text{ and } && 0 = X_1(0) = X_2(0) = X_3(0) = \ldots ,
\end{align*}
so let's go ahead and do that that. Plugging our Ansatz into
$$X'(x) = X(x)M(x)$$
and assuming it makes sense to differentiate term-by-term, we then get
$$ X_1'(x) + X_2'(x) + X_3'(x) \ldots = M(x) + X_1(x)M(x) + X_2(X)M(x) + \ldots.$$
A rather cheap way to make sure this equality holds is to match up terms on the left and right, i.e. require
\begin{align*}
X_1'(x) =  M(x) && X_2'(x) = X_1(x)M(x) && X_3'(x) = X_2(x)M(x) && \ldots
\end{align*}
This gives us a not particularly difficult recurrence relation to solve. Integrating (remember we put $X_n(0)=0$ for $n>0$), we get
\begin{align*}
X_1(x) &= \int_0^x M(x_1) \ dx_1 \\
X_2(x) &= \int_0^x X_1(x_2)M(x_2) \ dx_2 = \int_0^x \int_0^{x_2} M(x_1) M(x_2) \ dx_1 \ dx_2 
\end{align*}
and so on. The general term is:
$$ X_n(x) = \int\cdots\int_{0 \leq x_1 \leq \ldots \leq x_n \leq x} M(x_1) M(x_2) \cdots M(x_n) \ d x_1 \cdots d x_n.$$
Well OK, that's a bit imprecise. What does this mean when $x <0$? This is not a very serious issue, however, and can be completely eliminated by a simple change of variables.
Once we have the right formula, it's not all that difficult to do the analysis necessary to show that $X(x):=\sum_{n=0}^\infty X(x)$ indeed solves the IVP. One piece of relevant information is that the Euclidean volume of the $n$-simplex  $0 \leq x_1 \ldots \leq x_n \leq x$ is $x^n/n!$ as follows from noting that this simplex is, up to measure zero, a fundamental domain for the action of the symmetric group of order $n$ on the cube $[0,x]^n$  which permutes coordinates. This gives you estimates like
$$\|X_n(x)\| \leq \frac{|x|^n}{n!} \sup_{0 \leq y \leq x} \|M(y)\|^n$$
which allow you, for $x$ confined to any fixed bounded interval, to dominate the series $\sum_{n=0}^\infty X_n(x)$ by an exponential series.
