In order to solve a system of nonlinear rate equations one can try to use some soliton methods. The procedure involves three steps: (1) Find a `Lax representation' where all the kinetic variables are combined into a single matrix $\rho$, all the kinetic constants are encoded in a matrix $H$; (2) find a Darboux-Backund dressing transformation for the Lax representation $i\dot \rho=[H,f(\rho)]$, where $f$ models a time-dependent environment; (3) find a class of seed solutions $\rho=\rho$ that lead, via a nontrivial chain of dressings $\rho\to \rho\to \rho\to\dots$ to new solutions, difficult to find by other methods. The latter step is not a trivial one since a non-soliton method has to be employed to find an appropriate initial $\rho$. Procedures that lead to a correct $\rho$ have been discussed in the literature only for a limited class of $H$ and $f$. Here, we develop a formalism that works for practically any $H$, and any explicitly time-dependent $f$. As a result, we are able to find exact solutions to a system of equations describing an arbitrary number of species interacting through (auto)catalytic feedbacks, with general time dependent parameters characterizing the nonlinearity. Explicit examples involve up to 42 interacting species.