I have a nonlinear system of ODEs with known constant coefficients $A, B, C, D, E, F, M$:
\begin{align}
&\dot{n}(t)=-An(t)+Bm(t)n(t)+Cm(t) \\
&\dot{m}(t)=-Bm(t)n(t) + (M-m(t))R_0 \sin{\omega t} +Dm(t)+En(t)+F,
\end{align}
with driving amplitude $R_0$. I know from numerical simulations that solutions of $n(t)$ and $m(t)$ are periodic with the same period $T=2\pi/\omega$, which has the same period as my driving, i.e., $n(t)=n(t+T)$ and $m(t)=m(t+T)$. The problem is that these solutions are highly nonlinear, and only for small driving amplitude $R_0$ I can approximate them with $\cos(\omega t)$ term. I have figured out the solution, if I assume that $n$ and $m$ oscillate with a function $\cos(\omega t + \phi)$, but that is valid only when $R_0$ is small. I would like to somewhat obtain solution when $R_0$ gets larger, because then the solution does not oscillate like a $\cos (\omega t)$ anymore, but something like this:
Blue is the solution of $n(t)$, red is solution of $m(t)$, x-axis is time.
Is there a way to find a solution for general values of $R_0$ by utilizing the periodicity condition?