$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$ $\newcommand{\sgn}{\operatorname{sgn}}$ $\newcommand{\rank}{\operatorname{rank}}$
For $a > 0$ consider the nonlinear problem \begin{align} &u''+u-\varepsilon (u'-\frac{1}{3}u'^3)=0, \label{eq-8.4.1}\\ &u(0,\varepsilon)=0, \qquad u'(0,\varepsilon)=2a. \label{eq-8.4.2} \end{align} To remove secularities at leading-order, it is suffcient to introduce a single slow variable, in which case the two-variable expansion and the derivative expansion methods (with $N = 1$) are equivalent. So, $\tau=\varepsilon t$, \begin{equation} u(t,\varepsilon)\sim w(t,\tau,\varepsilon)=w_0(t,\tau)+\varepsilon w_1 (t,\tau). \label{eq-8.4.3} \end{equation} Then \begin{align} &\varepsilon^0:\ \left\{\begin{aligned} &\frac{\partial ^2w_0}{\partial t^2}+w_0=0,\\ &w_0(0,0)=0,\\ &\frac{\partial w_0}{\partial t}(0,0)=2a; \end{aligned}\right. \label{eq-8.4.4}\\[4pt] &\varepsilon^1:\ \left\{\begin{aligned} &\frac{\partial ^2w_1}{\partial t^2}+w_1= -2\frac{\partial ^2w_0}{\partial t\partial\tau} +\frac{\partial w_0}{\partial t}- \frac{1}{3}\bigl(\frac{\partial w_0}{\partial t}\bigr)^3,\\ &w_1(0,0)=0,\\ &\frac{\partial w_1}{\partial t}(0,0)=0. \end{aligned}\right. \label{eq-8.4.5} \end{align} The solution to problem (\ref{eq-8.4.4}) could be rewritten as \begin{align} &w_0(t,\tau) = A(\tau)\sin\bigl(t+\theta(\tau)\bigr),\label{eq-8.4.6}\\ &A(0)=2a>0,\quad\theta(0)=0.\label{eq-8.4.7} \end{align} Then equation in (\ref{eq-8.4.5}) is \begin{multline} \frac{\partial ^2w_1}{\partial t^2}+w_1=\\ \bigl(A-2A'-\frac{1}{4}A^3\bigr) \cos (t+\theta) +2A\theta' \sin(t+\theta) -\frac{1}{12} A^3\cos (3(t+\theta)) \label{eq-8.4.8} \end{multline} where we have expressed the right-hand side directly in terms of Fourier harmonics using the relation $\cos^3(t)=\frac{1}{4}\cos(3t)+\frac{3}{4}\cos(t)$.
To avoid secularities we equalize coefficients here to $0$: \begin{align} &A-2A'-\frac{1}{4}A^3=0,\label{eq-8.4.9}\\ &2A\theta'=0.\label{eq-8.4.10} \end{align} We can integrate (\ref{eq-8.4.9}): $8A'=A(4-A^2)$ and \begin{gather*} \int d\tau= \int \frac{-8\,dA}{A(A^2-4)}= \int \Bigl(\frac{2}{A}-\frac{2A}{A^2-4} \Bigr)\,dA= \ln \Bigl(\frac{A^2}{A^2-4}\Bigr)+\ln \alpha. \end{gather*} Then $(A^2-4)/A^2=\alpha e^{-\tau}$ with $\alpha=(a^2-1)/a^2$ because $A(0)=2a$. Therefore \begin{equation} A(\tau)=\frac{2a}{\sqrt{a^2-(a^2-1)e^{-\tau}}}>0. \label{eq-8.4.11} \end{equation} Then \begin{align} &u(t,\varepsilon)\sim\frac{2a\sin(t)}{\sqrt{a^2-(a^2-1)e^{-\varepsilon t}}},\label{eq-8.4.12}\\ &u'(t,\varepsilon)\sim\frac{2a\cos(t)}{\sqrt{a^2-(a^2-1)e^{-\varepsilon t}}}. \label{eq-8.4.13} \end{align} because $\theta(\tau)$ is constant, equal $\theta(0)=0$.
Since $A(\tau)\to 2$ as $\tau\to +\infty$ we see that the asymptotic solution approaches a limit cycle as $t\to +\infty$.
Remark 1. However, this is true for asymptotic solution $w_0 (t,\varepsilon t) )+ \varepsilon w_1 (t,\varepsilon t)$ which is approximates exact solution $u_\varepsilon (t)$ only as $\varepsilon^2 t \ll 1$ with an error $O(\varepsilon ^2 t)$. Using properties of 2-dimensional dynamical systems one can prove that exact has a limit cycle which is contained in $C\varepsilon$-vicinity of this circle.
To derive more exact result consider (\ref{eq-8.4.8}) which becomes
\begin{gather*}
\frac{\partial ^2w_1}{\partial t^2}+w_1=
-\frac{1}{12}A^3(\tau) \cos (3t)
\end{gather*}
which with initial conditions implies that
\begin{gather*}
w_1 = \frac{1}{96}A^3 (\tau) \cos (3t).
\end{gather*}
However, we need also to replace (\ref{eq-8.4.3}) by
\begin{equation}
u(t,\varepsilon)\sim w(T,\tau,\varepsilon)=w_0(T,\tau)+\varepsilon w_1 (T,\tau)
\label{eq-8.4.14}
\end{equation}
with $T=t +\nu \varepsilon^2 t$. Then all previous analysis remains true with $t$ replaced by $T$ and also we have an equation
\begin{gather*}
\frac{\partial^2 w_2}{\partial T^2} =
-\nu \frac{\partial^2 w_0}{\partial T^2} +
\Bigl(1-\bigl(\frac{\partial w_0}{\partial T}\bigr)^2\Bigr)\frac{\partial w_1}{\partial T}.
\end{gather*}
Plugging $w_0,w_1$ and selectiong $\nu$ to get rid of secular term we can construct $w_2$ and thus
$u_\varepsilon (t)$ modulo $O(\varepsilon^3 t)$ as $\varepsilon^3 t\ll 1$. Also, without $\varepsilon^2 w_2$ the error is $O(\varepsilon^3t + \varepsilon^2)$.
Then the limit cycle of the exact solution is contained in $C\varepsilon^2$ vicinity of $(w_0+\varepsilon w_1, w'_0+\varepsilon w'_1)$ with $T$ replaced by $t$ and $A=2$.
Similar construction can be applied to Van der Pol oscillator \begin{align} &u''+u -\varepsilon (1-u^2)u' =0 ,\label{eq-8.4.15}\\ &u(0)=1,\qquad u'(0)=0; \label{eq-8.4.16} \end{align} and while in (\ref{eq-8.4.5}) the last term in the right hand expression is replaced by $-w_0^2 \frac{\partial w_0}{\partial t}$ and in (\ref{eq-8.4.8}) the coefficient in the last term differs, where we used $\sin^2(t)\cos(t)=\frac{1}{4}\cos(t)-\frac{1}{4}\cos(3t)$, (\ref{eq-8.4.9})--(\ref{eq-8.4.14}) remain unchanged.
Remark 2.