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.
The solution to the initial value problem
$$
\dot{x}=-x+x^2,\quad x(0)=x_0
$$
is
$$
x(t)=\frac1{1-(1-1/{x_0})e^t}.
$$
Its derivative is
$$
\dot x(t)= \frac{(1-1/x_0)e^t}{(1-(1-1/{x_0})e^t)^2}.
$$
Note that the denominator of this fraction can't be equal to zero if $t\geq 0$.
Let $x_0\in(-1,1)\setminus\{0\}$. We have $1-\frac1{x_0}<0$ if $x_0\in(0,1)$ and $1-\frac1{x_0}>0$ if $x_0\in(-1,0)$, thus $\dot x(t)$ is positive ($x(t)$ is monotonically increasing) if $x_0\in(-1,0)$ and $\dot x(t)$ is negative ($x(t)$ is monotonically decreasing) if $x_0\in(0,1)$. Since the solution to the initial value problem cannot cross the equilibrium point $x=0$, its norm $\|x(t)\|$ decreases monotonically for all $t\ge 0$ and $x_0\in(-1,1)\setminus\{0\}$.
This means that we can choose $\delta(\varepsilon)= \min(\varepsilon,1)$ in the stability definition.
Best Answer
Multiplying the second equation by $y$ and using the first equation we have $$y\dot{y}=-2xy-p(t)y^4=-x\dot{x}-p(t)y^4$$ i.e. as correctly pointed in the comments by @Evgeny $$\frac{d}{dt}[x^2+y^2]=-2p(t)y^4\leq -4 y^4$$
Thus $x,y$ are bounded and the Lyapunov-like function $V=x^2+y^2$ is decreasing and lower bounded (from zero). Hence, $V$ converges to some constant $V_{\infty}\geq 0$. If we integrate the above inequality over $[0,\infty)$ we obtain $$V_{\infty}-V(0)\leq -4\int_0^{\infty}{y^4(s)ds}$$ which yields $$\int_0^{\infty}{y^4(s)ds}\leq \frac{1}{4}V(0)$$
If we assume a bounded $p(t)$ then the boundedness of $x,y$ and the state equations result in the boundedness of $\dot{x},\dot{y}$.
A continuous differentiable function with bounded derivative is uniformly continuous.
You can use now Barbalat's lemma to prove convergence. Barbalat's lemma states that if
i) $\int_0^{\infty}{\phi(t)dt}$ exists and is finite and
ii) $\phi(\cdot)$ is uniformly continuous
then $\lim_{t\rightarrow\infty}\phi(t)= 0$.
From Barbalat lemma we then have that $\lim_{t\rightarrow\infty}y(t)=0$ and since $V$ converges there also exists some $x^*$ (with $V_{\infty}={x^*}^2$) such that $\lim_{t\rightarrow\infty}x(t)=x^*$.
Edit to remove the upper bounded $p(t)$ restriction: I will prove now a variation of Barbalat lemma that does not need uniform continuity of $\phi(\cdot)$ but only an upper bounded derivative.
Barbalat lemma variation: Let $\phi:\mathbb{R}\rightarrow\mathbb{R}_+$ continuous differentiable nonnegative function. If
i) $\int_0^{\infty}{\phi(t)dt}$ exists and is finite and
ii) $\dot{\phi}(\cdot)$ is upper bounded
then $\lim_{t\rightarrow\infty}\phi(t)= 0$.
Proof: The proof uses a contradiction argument. Assume the opposite, then there exists a constant $k_1>0$ and a sequence of times $\{T_i\}$ with $\lim_{i\rightarrow\infty}T_i=\infty$ such that $$\phi(T_i)\geq k_1$$ For $t\leq T_i$ we have from the mean value theorem that $$\phi(t)=\phi(T_i)-\dot{\phi}(\theta t +(1-\theta)T_i)(T_i-t)$$ for some $\theta\in(0,1)$. Since $\dot{\phi}$ is upper bounded there exists some $c$ such that $\dot{\phi}(t)\leq c$ for all $t\geq 0$ and $$\phi(t)\geq \phi(T_i)-c(T_i-t)\qquad \forall t\leq T_i$$ From the above relationship we have that $$\phi(t)\geq \frac{k_1}{2}\qquad \forall t\in\left[T_i-\frac{k_1}{2|c|},T_i\right]$$ Therefore $$\int_{T_i-\frac{k_1}{2|c|}}^{T_i}{\phi(t)dt}\geq \frac{k_1^2}{4|c|} $$ Thus $\int_0^{t}{\phi(\tau)d\tau}$ cannot converge to a finite limit as $t\rightarrow \infty$ which is the desired contradiction.
In our example $\phi(t):=y^4$ and its derivative is upper bounded since $$\dot{\phi}(t)=4y^3\dot{y}=-8xy^3-4p(t)y^6\leq -8xy^3$$ and $x,y$ are bounded. The proposed variation of Barbalat lemma yields now the desired $\lim_{t\rightarrow\infty}y(t)=0$. This completes the proof and indeed there is no need for an upper bound of $p(t)$.