One can multiply the equation with $2x'$ and integrate to get a first integral that has bounded level sets around the origin. This gives raise to closed orbits and periodic solutions.
It is expected that the period will also depend on $\epsilon$, therefore one tries the solution form $x(t)=C(\omega(\epsilon)t,\epsilon)$. The first argument of $C$ or the product inside gets the name $\theta$, and by abuse of notation $C$ is replaced with $x$, even though it is now a different function. This way one gets
$$
x''(s)=\omega(\epsilon)^2\frac{\partial^2 C}{\partial\theta^2}.
$$
It is $C(\theta,\epsilon)$ that gets expanded as $x_0(\theta)+\epsilon x_1(\theta)+...$ which leads to the mentioned detailed equation
\begin{multline}
(1+ϵω_1+ϵ^2ω_2+...)^2(x_0''(θ)+ϵ x_1''(θ)+ϵ^2 x_2''(θ)+...)\\+(x_0(θ)+ϵ x_1(θ)+ϵ^2 x_2(θ)+...)\\=ϵ(x_0(θ)+ϵ x_1(θ)+ϵ^2 x_2(θ)+...)^2
\end{multline}
and by comparing coefficients to the system
\begin{align}
x_0''(θ)+x_0(θ)&=0\\
2ω_1x_0''+(x_1''+x_1)&=x_0^2\\
(2ω_2+ω_1^2)x_0''+ω_1x_1''+(x_2''+x_2)&=2x_0x_1
\end{align}
etc. So one finds $x_0(θ)=\cos(θ)$.
The aim is to fix the parameters and component functions to avoid non-periodic resonance terms as one would get from $2ω_1x_0$ in the second equation. So set $ω_1=0$. Thus remains
$$
x_1''(θ)+x_1(θ)=x_0(θ)^2=\frac12(1+\cos(2θ))
$$
From the method of undetermined coefficients the general solution form is known as
$$
x_1(θ)=A+B\cos(2θ)+C\cos(θ)+D\sin(θ)
$$
Now one can reduce the system by directly satisfying the homogeneous initial conditions $x_1(0)=x_1'(0)=0$ by shifting some amount from the complementary solution to the particular solution
$$
x_1(θ)=A(1-\cos(θ))+B(\cos(2θ)-\cos(θ))
$$
Inserting into the DE results in $A=\frac12$, $-3B=\frac12$.
In the third equation one can replace $x_0''$ with $-x_0$ to get
\begin{align}
x_2''+x_2&=2ω_2x_0+2x_0x_1\\&=(2ω_2+1-\tfrac16)\cos(θ)-\tfrac16(1+\cos(2θ))-\tfrac1{6}\cos(3θ)
\end{align}
Thus $ω_2=-\frac5{12}$ to avoid non-periodic resonance terms and so on.
Trying out how this looks like numerically, integrate with $ϵ=0.05$ as example.
This looks like $x=y+ϵu$, $u(0)=0$, that is, in your expansion the $X_0$ component does not oscillate, the oscillation happens in $X_1$, so what you computed is actually correct, you just need to advance to the next step.
Alternatively transform first the system without series expansion. Eliminate $x$ using $y-x = -ϵu$ and $x-2y = -(y-ϵu)$, and in the derivative
$$
ϵ v \overset{Df.} = \dot x =\underbrace{-ϵy(y-ϵu)}_{=\dot y} + ϵ\dot u
\\
\implies
\dot u = v + y(y-ϵu).
$$
Now insert into the first equation
$$
\dot v = \frac{1}{ϵ}\ddot x = \frac{1}{ϵ}(y-x)=-u ,~~v(0)=0
$$
and as already used, the second equation transforms to
$$
\dot y = - ϵy(y-ϵu)
$$
Plotting this system shows all components at the same scale
One could now switch to a systematic perturbation expansion, but let's go just one step further in the manual transformation. Set $w=v+y^2$, $w(0)=1$, then
\begin{align}
\dot u &= w - ϵyu\\
\dot w &= -u - 2ϵy^2(y-ϵu)\\
\dot y &= - ϵy(y-ϵu)
\end{align}
Thus with $z=u+iw$
$$
\dot z=-iz-ϵy(Re(z)+i2y(y-ϵRe(z))),\\ z(0)=i\implies z=ie^{-it}+O(ϵ)
$$
At level 0 this properly decouples as $u=\sin t$, $w=\cos t$, $y=1$. For a two-time-scales approach you would next introduce the slow variable $T=ϵt$ and equip the solution with slowly changing constants
$$
z_0(t,T)=A(T)e^{-it}+B(T)e^{it}, y_0(t,T)=C(T) \\
A(0)=i,~~B(0)=0,~~ C(0)=1.
$$
On the next level, with $z=z_0+ϵz_1+...$, etc.
\begin{alignat}{1}
\dot z_1 &= -iz_1 - C(T)Re(A(T)e^{-it}+B(T)e^{it})-i2C(T)^3&-A'(T)e^{-it}-B'(T)e^{-it}\\
\dot y_1 &= - C(T)^2 -C'(T)
\end{alignat}
The secular terms in the first equation are those that contain $e^{-it}$. To avoid resonances, set them to zero
$$
0=\frac{C(T)}2(A(T)+\overline{B(T)}+A'(T)
$$
The remaining equation gives
$$
\dot z_1+iz_1 = + \overline{A'(T)}e^{it} - 2iC(T)^3 - B'(T)e^{it}
$$
It appears natural to demand $\overline{A'(T)} - B'(T) = 0$, giving $A(T)-\overline{B(T)}=i$ and using $C=-C'/C$ leading to
$$
A(T)=\frac i2(1+C(T)), \\ B(T)=\frac i2(1-C(T)).
$$
It remains to solve
$$
\dot z_1+iz_1=-2iC^3, ~~ z_1(0)=0,\\
z_1 = -2C^3 + 2C_1(T) e^{-it}, ~~~C_1(0)=1.
$$
So in reverse the steps are
i = 1j
T = eps*t
C = 1/(1+T); A = 0.5*i*(1+C); B = 0.5*i*(1-C)
z_0 = A*exp(-i*t) + B*exp(i*t) # = -C*sin(t) + i*cos(t)
z_1 = -2*C**3*(1-exp(-i*t))
u = z_0.real + eps*z_1.real # = -C*sin(t) + 2*eps*(cos(t)-C**3)
w = z_0.imag + eps*z_1.imag # = +cos(t) -2*eps*sin(t)
y = C
v = w - y**2
dx = eps*v
x = y + eps*u
giving the plot
This looks qualitatively correct, quantitatively only in the short term.
Best Answer
You also need to insert the expansion of $Ω(ϵ)=1+ϵΩ_1+ϵ^2Ω_2+...$ so that you get $$ (1+2ϵΩ_1+..)(X_0''+ϵX_1''+..)+(X_0+ϵX_1+...)=ϵ(X_0^2+2ϵX_0X_1+...) $$ Comparing coefficients then gives in the first two orders of $ϵ$ $$ X_0''+X_0=0\\ 2Ω_1X_0''+X_1''+X_1=X_0^2\\ $$ As $X_0^2$ will have frequencies $0$ and $2$, there are no secular terms to compensate, so set $Ω_1=0$. Then in the $ϵ^2$ terms you get $$ 2Ω_2X_0''+X_2''+X_2=2X_0X_1\\ $$ The product $X_0X_1$ consist of terms with frequencies $2\pm1$ and $0+1$, so there will be secular terms of frequency $1$ that can be compensated by the $X_0''$ term, giving $Ω_2$ a non-zero value.
In detail, $X_0=\cos(θ)$, resulting in $X_0^2=\frac12(1+\cos(2θ)$, which gives per "undetermined coefficients" $$X_1=A(1-\cos(θ))+B(\cos(2θ)-\cos(θ)),$$ with $A=\frac12$ and $-3B=\frac12$. Now $$ 2X_0X_1=\frac12(2\cosθ-1-\cos(2θ))-\frac16(\cos(3θ)+\cosθ-1-\cos(2θ)) $$ thus $-2Ω_2=1-\frac16=\frac56$.