Assume $y_1<y_2$ and $f(x,y_1)>0$, $\ f(x,y_2)<0$ for all $x$. By general principles about ODE's any solution $x\mapsto \phi_\eta(x)$ starting at a point $(0,\eta)$, $\ y_1\leq\eta\leq y_2$, will finally leave the rectangle $R:=[0,w]\times[y_1,y_2]$. But it cannot do so along the horizontal edges of this rectangle. It follows that the solution $\phi_\eta$ will pass through a point $(w,\eta')$, $\ y_1\leq\eta'\leq y_2$, on the right edge of $R$. In this way a so-called Poincaré map
$$\Phi:\quad [y_1,y_2]\to[y_1,y_2],\qquad \eta\mapsto \eta'=:\Phi(\eta)$$
is defined. Again by general principles this $\Phi$ is continuous. By Brouwer's fixed point theorem (or using the intermediate value theorem) it follows that $\Phi$ has a fixed point $\eta_*\in[y_1,y_2]$. The solution $\phi_{\eta_*}$ starting at $(0,\eta_*)$ is then periodic.
If instead of the initial assumption on $f$ we have $f(x,y_1)<0$, $\ f(x,y_2)>0$ we start the argument at $x=w$ and proceed to the left.
In the example $y'+p(x) y=q(x)$ we have
$$f(x,y)=p(x)\left({q(x)\over p(x)} -y\right)\ .$$
As $p$ and $q$ are periodic and continuous, and $p(x)\ne0$ for all $x$ there is an $M>0$ such that
$$-M<{q(x)\over p(x)}<M\qquad\forall x\ .$$
It follows that by choosing $y_1:=-M$, $\ y_2:=M$ we can fulfill the assumptions of the "theorem".
This is not a satisfactory proof of the OP's question. It's a sketch of a proof of a slightly weaker claim (see assumptions below). However, since this question hasn't gotten attention from an expert, as I hoped for in a comment above, then I'll write this down as the best that I can do at the moment. There may be counterexamples to my argument, but I think they would have to be pathological.
One important point: I'm interpreting that the question requires that every point be in a strictly periodic orbit, in the sense that for every initial condition $x(0)$, there will be some $t > 0$ such that $x(t) = x(0)$. I think that this is the sense in which the OP meant that "all solutions...are periodic." (A physicist like me might informally say that all solutions are periodic if they asymptotically approach one or more periodic attractors, which is a completely different situation.)
I'll make the following assumptions that are stronger than the OP's question. These can probably be relaxed, but I have not spent any time trying to do so.
- $f$ is differentiable (not just Lipschitz).
- The Jacobian of $f$ at $x=0$, which I'll call $J$, is not identically zero.
Proof: I will try to prove that your professor's picture is essentially correct, i.e., that initial conditions near the origin follow orbits that remain near the origin. Sufficiently close to the origin,
$$ \dot{x} = f(0) + J x + \dots \approx J x $$
since $f(0)=0$. This linearization has the solution
$$ x(t) = e^{Jt} x(0) $$
as long as $x(t)$ is sufficiently close to $0$.
The eigenvalues of $J$ are the local Lyapunov exponents. Note that these can be complex, but since $J$ is real, they must come in complex conjugate pairs; linear combinations of these are used to create real solutions for $x(t)$. (I will ignore this subtlety below but it doesn't change anything.) In the next paragraph, I try to show that $\text{Re}(\lambda) = 0$ for all eigenvalues $\lambda$.
Suppose that there is some eigenvector $v$ whose associated eigenvalue $\lambda$ has a negative real component, i.e., $\text{Re}(\lambda) < 0$. Suppose we let $x(0) = v$. Then $x(t) = v e^{\lambda t}$, and $x(t)$ will asymptotically approach $0$, a fixed point. Specifically, $x(t)$ will never return to $x(0)$, which we assumed that it would. Note that the linearization of the original ODE is still valid within this entire neighborhood, so we can't say that somewhere along the path, it will get swept off in a different direction and complete a periodic orbit without going to $0$.
Dealing with the case that $\text{Re}(\lambda) > 0$ is slightly harder. This is because, after some time, $x(t)$ will leave the neighborhood where the linearization is valid, and after that we cannot say what will happen. But here I will employ a trick: Make time go backwards. Note that a periodic orbit must remain periodic whether time runs forward or backward. Now consider the original ODE,
$$ \dot{x} = f(x) $$
under the substitution $t \to -t$ (running time backwards). We get the ODE
$$ \dot{x} = -f(x) $$
which corresponds to taking $J \to -J$ and $\lambda \to -\lambda$ in the discussion above. So if we have $\text{Re}(\lambda) > 0$ when time is running forward, we can run time backward and apply the same argument as we did for $\text{Re}(\lambda) < 0$.
Note that because we can choose the initial condition to be exactly on the eigenvector, the above discussion applies even if $J$ has both positive and negative eigenvalues, because we are singling out the effect of a single eigenvalue.
We have shown that the cases where $x(t)$ move rapidly toward or away from the origin are excluded, so $x(t)$ must remain at approximately the same distance from the origin. In particular, the eigenvalues of $J$ must be pure imaginary, so $x(t)$ will be a sum of oscillating components. In particular, if $x(0) = \sum_i c_i v_i$ for eigenvectors $v_i$, then $x(t) = \sum c_i e^{\lambda_i t} v_i$. Obviously $|x(t)|$ can change, but it is bounded by the triangle inequality.
The main place where this argument needs more rigor is that one would need to do an $\epsilon$-$\delta$ analysis (or invoke relevant theorems) to validate that we really can ignore the effect of the higher-order terms in the Taylor series. But I think the basic outline is sound.
Best Answer
You found that along a solution the expression $$ \dot x^2+x^2(1+\tfrac12x^2)=R^2 $$ is constant. The left side can be seen a sum of two squares, so that the whole is a circle equation for radius $R$. Parametrize with the angle as $$\dot x=R\cos(u), ~~ x\sqrt{1+\frac12 x^2}=f(x)=R\sin(u).$$ Then taking the derivative of the second expression and comparing with the first gives $$ f'(x(t))\dot x(t)=R\cos(u(t))\dot u(t)\\ \implies f'(x(t))=\dot u(t) $$ for any motion where $\dot x$ is not constant zero.
Now $f(x)$ is a strictly increasing function, $f'(x)\ge 1$ is always positive. This means that $u(t)$ is also strictly increasing, and the phase curve $(x(t),\dot x(t))$ moves along the image of a circular motion, thus is periodic.
$$ f'(x)=\sqrt{1+\tfrac12 x^2}+\frac{x^2}{2\sqrt{1+\tfrac12 x^2}}=\frac{1+x^2}{\sqrt{1+\tfrac12 x^2}}\ge 1 $$ and $y=f(x)$ can be solved observing $\operatorname{sgn}(y)=\operatorname{sgn}(x)$ as follows $$ y^2=x^2+\tfrac12x^4\\ 1+2y^2=(1+x^2)^2\\ x^2=\sqrt{1+2y^2}-1=\frac{2y^2}{1+\sqrt{1+2y^2}}\\ \implies x=\frac{\sqrt2 y}{\sqrt{1+\sqrt{1+2y^2}}} $$