I have encountered a system of nonlinear differential equations that look and behave similar to the predator-prey model:
$$\begin {align*} \frac{dx}{dt} &= (1+x) \left( \frac{1+x}{y} - y \right) \\ \frac{dy}{dt} &= (1+y) \left( x - \frac{1+y}{x} \right) \end{align*} $$
Superficially, it resembles the predator-prey model. But also, the solution trajectories behave in much the same way. If you start with positive $x$ and $y$, then the solution curves are closed curves in the first quadrant.
I'm wondering if anybody knows any techniques to explicitly integrate these equations and find the functions $x(t)$ and $y(t)$ for arbitrary initial conditions (let's say $x(0) >0$ and $y(0) > 0$).
Extra Background: I happen to know already that that these equations are a Hamiltonian flow with respect to the Poisson bracket $\{x,y\} = xy$ for the Hamiltonian $$ H = x + y + \frac{1+y}{x} + \frac{1+x}{y} + \frac{1+x+y}{xy} $$ Indeed, that is how I came up with the differential equations above.
The solution $(x,y)$ seems to be periodic in time $t$. The harmonic balance method is a semi-analytic procedure that consists in looking for the Fourier coefficients of $x$ and $y$. Consider for instance the first-order truncated Fourier series \begin{aligned} x(t) &= a_0(x) + a_1(x)\cos(\omega t) + b_1(x)\sin(\omega t) \\ y(t) &= a_0(y) + a_1(y)\cos(\omega t) + b_1(y)\sin(\omega t) . \end{aligned} Then, inject this Ansatz in the system and compute the scalar products of each equation with $\cos(n\omega t)$, $\sin(n\omega t)$. This will give you the relationship between the coefficients of this truncated series. The process can be repeated up to an arbitrary number of harmonics, allowing the identification of the solution's Fourier series.