@DanieleTampieri actually, since the small-$t$ behavior is reasonably well-approximated by a single Runge-Kutta step (and numerical solvers remain stable up to $t \sim O(\lambda^{-1}+\mu^{-1})$), I’d be even more interested in finding approximately asymptotic behavior as $t$ gets large. I know $a\to 0$ and empirically I can see a few other things from the numerics but it’d be much better to have bounds or good approximations…

@DanieleTampieri thank you for that; yes any kind of insight that I can get into these would be helpful, including a function series expansion for small $t$. Do you have any thoughts or suggestions as to how I might find a suitable basis? Something like orthogonal polynomials? Would the orthogonality be defined over the entire domain $[0,\infty)$?

@CarloBeenakker thanks, yes I suspected that as a nonlinear system the default would be that it’s unsolvable, but the existence of a solution for the special case $x=y=0$ led me to wonder if anything could be done in the more general case