I found the following problem during my research (in cosmology):
given the functions $P_1(x),\ P_2(x),\ P_3(x)$, how to solve the system
\begin{equation} \frac{d\ P_m}{dx} = P^t \cdot A_m \cdot P - P_m\ (P^t \cdot B \cdot P)\,, \qquad \text{(1)} \end{equation}
where $P^t= (P_1,P_2,P_3)$ is the transposed vector of $P$, $B(x)$ and $A_1(x),\ A_2(x),\ A_3(x)$ are arbitrary (but known) matrices, and a "$\cdot$" means product between the vectors and the matrices. Here $P_m$ is the $m$-component of the vector $P$.
PERTURBATIVE SOLUTION: I can solve the system perturbatively. For instance, if I suppose that the matrices $A_m$ and $B$ are small. Then, given the initial condition $P_m(x=0)= q_m$, I can to first order set $P_m = q_m$ on the r.h.s of equation (1) to obtain:
\begin{equation} P_m = q_m + q^t \cdot \left(\int dx\ A_m \right) \cdot q - q_m\ \left[q^t \cdot \left( \int dx\ B \right) \cdot q\right]\,, \quad \text{(2)} \end{equation}
and use this solution back into (1) to obtain a second order solution, and so on.
Is there a way to solve the system without going into perturbation theory?