Finding $T$-period solutions of nonlinear first order ODEs Problem:
There is a set of coupled ODEs, where the time derivatives are equal to quartic polynomials. Here is an example:
\begin{equation}
\begin{aligned}
\frac{dx_1}{dt} &= x_2 x_3 (y_1^2 + y_3^2 + 1) \\
\frac{dx_2}{dt} &= -x_1 x_3 (y_2^2 + y_3^2 + 1) \\
\frac{dx_3}{dt} &= x_1 x_2 (y_1^2 - y_2^2) \\
\end{aligned}
\qquad
\begin{aligned}
\frac{dy_1}{dt} &= y_2 y_3 (x_1^2 + x_3^2 + 1)  \\
\frac{dy_2}{dt} &= -y_1 y_3 (x_2^2 + x_3^2 + 1) \\
\frac{dy_3}{dt} &= y_1 y_2 (x_1^2 - x_2^2) \\
\end{aligned}
\end{equation}
Periodic (numerical) solutions with period $T$ are wanted, i.e. $\vec{x}(0) = \vec{x}(T)$, $\vec{y}(0) = \vec{y}(T)$. The initial condition is not known, otherwise one would have easily solved it numerically via Rugge-Kutta. However, the boundary condition can be vaguely posed as follows. We start from the linear problem
\begin{equation}
\begin{aligned}
\frac{dx_1}{dt} &= x_2 x_3  \\
\frac{dx_2}{dt} &= -x_1 x_3 \\
\frac{dx_3}{dt} &= 0 \\
\end{aligned}
\qquad
\begin{aligned}
\frac{dy_1}{dt} &= y_2 y_3   \\
\frac{dy_2}{dt} &= -y_1 y_3 \\
\frac{dy_3}{dt} &=  \\
\end{aligned}
\end{equation}
with constant $\vec{x}(t) = (x_1, 0, 0)$ and $\vec{y}(t) = (0, 1, 0)$. Then gradually turn on the non-linear terms. There may be an iterative process leading to the periodic solutions of the non-linear equations.
Analytic Attempt:
For the example problem, $\vec{y} = (0, 1, 0)$ solves the $y$ equations without knowing $x$. Then the $\vec{x}$ equation reduces to
\begin{equation}
\begin{aligned}
\frac{dx_1}{dt} &= x_2 x_3 \\
\frac{dx_2}{dt} &= -2x_1 x_3 \\
\frac{dx_3}{dt} &= -x_1 x_2
\end{aligned}
\end{equation}
Possible solutions for $T = \infty$ can be
\begin{equation}
\begin{aligned}
x_1 &= \tanh ( \sqrt{2} t ) \\
x_2 &= \frac{\sqrt{2}}{\cosh( \sqrt{2} t  )}\\
x_3 &= - \frac{1}{\cosh( \sqrt{2} t )}
\end{aligned}
\end{equation}
             
             
             
             

In this solution $x_1$ goes up from $0$ to $1$, I imagine in a periodic solution, it goes down symmetrically to $0$ on the other end.
Numerically calculation suggests that perturbing the initial conditions can produce periodic solutions, but how the initial condition controls the period is not clear.
Question:
This is probably an analytic solution for the example problem. However, the actual sets of ODEs contains $3n$ variables, so I seek to find a numerical approach.
Is there a variation method to find periodic solutions with unspecified boundary conditions?
Suggestions of other methods are also very welcome.
Thank you!

Just a suggestion: there are two conserved quantities, $$x_1^2+x_2^2-x_3^2$$ and similarly for the $y$s. That brings the problem from 6 dimensions down to 4. If you could find 2 more, there would be a lot of theory available to help.
