Solution (source code)

= Solution

Treat all transformation coefficients as $O(\lambda)$, so products of them with $\lambda$ are second order. Substitution into the <Hamiltonian> gives its harmonic part plus
$$
(m\omega^2\alpha_1+\lambda)X^3+
(2m\omega^2\alpha_2+\beta_1/m)X^2P+
(m\omega^2\alpha_3+2\beta_2/m)XP^2+
(\beta_3/m)P^3.
$$
Set these four coefficients to zero and impose the canonical constraints. The result is
$$
\boxed{\alpha_1=-\frac\lambda{m\omega^2},\quad
\alpha_2=0,\quad\alpha_3=-\frac{2\lambda}{m^3\omega^4},\quad
\beta_1=0,\quad\beta_2=\frac\lambda{m\omega^2},\quad\beta_3=0.}
$$
Thus the <first-order canonical removal of a cubic oscillator perturbation> is
$$
\boxed{x=X-\frac\lambda{m\omega^2}X^2-\frac{2\lambda}{m^3\omega^4}P^2,
\qquad p=P+\frac{2\lambda}{m\omega^2}XP,}
$$
up to $O(\lambda^2)$. It yields $K=P^2/(2m)+m\omega^2X^2/2+O(\lambda^2)$. The inverse map, if old <initial data> are specified, is $X=x+\lambda x^2/(m\omega^2)+2\lambda p^2/(m^3\omega^4)$ and $P=p-2\lambda xp/(m\omega^2)$ to the same order. There is no time-dependent generating term because the map is autonomous.