So far this is what I have:
First off, let $\dot x = y$ and $\dot y = -x-x^3$. Then the DE can be written as $$\frac{y^2}{2}+\frac{x^2}{2}+\frac{x^4}{4} = H(x,y)$$
This is Newtonian, so consider the following:
$$U(x) = \frac{x^2}{2} + \frac{x^4}{4}$$ $$U'(x) = x+x^3$$
This shows that $U'(x) = 0$ means that $x = 0$. Now
$$U''(x) = 1 + 3x^2$$
So $U''(0) = 1 > 0$ which implies that we have a min at $x = 0$. This implies that we have a center at $x =0$.
Is there more that I can add to this without explicitly using software? This is a question from a former PhD qualifying exam, so I'm trying to figure out what they may be expecting.
As additional information you can approximate the period of small oscillatins. Consider the solution curve through $(x,y)=(x_0,0)$, $x_0$ small. Then with the Hamiltionian we get $$ y^2+x^2+\frac12x^4=x_0^2+\frac12x_0^4 \\ \dot x^2=y^2=(x_0^2-x^2)\left(1+\frac12(x_0^2+x^2)\right). $$ One quarter of the period can be computed as the integral \begin{align} \frac{T}4&=\int_0^{x_0}\frac{dx}{\sqrt{x_0^2-x^2}\sqrt{1+\frac12(x_0^2+x^2)}}\\ &=\int_0^1\frac{dz}{\sqrt{1-z^2}\sqrt{1+\frac{x_0^2}2(1+z^2)}} \\[1em] &=\int_0^1\frac{dz}{\sqrt{1-z^2}\sqrt{1+x_0^2}\sqrt{1-\frac{x_0^2(1-z^2)}{2(1+x_0^2)}}} \\ &=\sum_{k=0}^\infty\binom{2k}{k}\frac{x_0^{2k}}{8^k(1+x_0^2)^{k+\frac12}}\int_0^1(1-z^2)^{k-\frac12}\,dz \\ &=\frac1{\sqrt{1+x_0^2}}\int_0^1\frac{dz}{\sqrt{1-z^2}} +\frac{x_0^2}{4(1+x_0^2)^{\frac32}}\int_0^1\sqrt{1-z^2}\,dz+O(x_0^4)\\[1em] &=\frac\pi{2\sqrt{1+x_0^2}}+\frac{\pi x_0^2}{16(1+x_0^2)^{\frac32}}+O(x_0^4). \end{align} The series expansion was done using the binomial series $$ (1-u)^{-1/2}=\sum\binom{-1/2}k(-u)^k=\sum\binom{2k}{k}\left(\frac u4\right)^k $$