Solution (source code)

= Solution

Insert the harmonic variables into the forward map and reduce the squared trigonometric functions. This gives
$$
\boxed{x(t)=A\cos\vartheta-
\frac{3\lambda A^2}{2m\omega^2}+
\frac{\lambda A^2}{2m\omega^2}\cos2\vartheta+O(\lambda^2),}
$$
$$
\boxed{p(t)=-m\omega A\sin\vartheta-
\frac{\lambda A^2}{\omega}\sin2\vartheta+O(\lambda^2).}
$$
Indeed $p=m\dot x$ at this order. The <second harmonic and mean displacement of a weakly cubic oscillator> show asymmetry of the motion: for positive $\lambda$ the mean shifts towards negative $x$. There is \b[no first-order <frequency> shift], although the waveform is no longer purely sinusoidal.

The small parameter is $|\lambda|A/(m\omega^2)\ll1$. A purely cubic perturbation is unbounded in one direction, so these formulas describe small local <oscillations> around the stable origin, not globally confined motion at arbitrary <energy>. Amplitude $A$ labels the transformed orbit and must be inferred using the inverse map if physical initial position and <momentum> are supplied.