Integrable nonlinear oscillators with polynomial invariants: construction, Poincare geometry, and an analytic stability boundary

Starting from the nonlinear ODE $z'' + f(t)\,z + g(t)\, z^{m}=0$ with $m>1$, we show that after a suitable normal-form reduction of any Hill equation one may, without loss of generality, fix the linear part as $f(t)\equiv \omega^{2}$ (with $\omega>0$ constant). For the class $z''+\omega^{2}z+g(t)\, z^{m}=0$ with $m>1$, our goal is to compile a catalogue of all possible integrable cases. We restrict attention to integrals that are polynomial in the variables $z$ and $p=z'$. The Hamiltonian does not provide such an integral because it is explicitly time dependent. Instead, we search for invariants that are quadratic in $p=z'$. We show that such invariants exist precisely when $\alpha_2(t):=g(t)^{-2/(m+3)}$ satisfies the linear third-order ODE $\alpha_2''' + 4\omega^2 \alpha_2'=0$. This yields the three-parameter solution $g(t)=[a_0+a_1\cos(2\omega t)+a_2\sin(2\omega t)]^{-(m+3)/2}$. For $m=2$ this reproduces the trigonometric structure with exponent $-5/2$ found in Hagel--Bouquet (1992). In addition we present a detailed stability analysis based on the invariant using Poincar\'e sections and find full agreement with numerical simulations.