Breaking the circle
Obviously, it is wrong to set $A=B=0$ as the method would prescribe. Meaning that you made an error in setting up the method, see below. Let's first look what happens when applying the method of undetermined coefficients
$$
x_1=Cτ\cosτ+Dτ\sinτ+\text{terms of double frequency}
$$
that is, you get resonance terms but no new contribution to the original homogeneous solution. Solving this gives $C=\frac{Bδ}2$, $D=-\frac{Aδ}2$ so that in total
$$
x_1=\frac{Bδ}2τ\cosτ-\frac{Aδ}2τ\sinτ+\frac{A}6(\cos2τ-3)+\frac{B}6\sin2τ
$$
On finite time spans these terms in $x_0+ϵx_1$ can be combined per Taylor to give an approximation using the sine and cosine of $(1+\frac{δϵ}2)τ$. This implies that the slow time should have been chosen as $T=δϵt$ or simply $T=ϵt$. Which does not lead toward the claim.
On solving the problem
You want to explore $a$ of the form $1+cϵ^2$. You can not do this by setting $δ=cϵ$, as all other quantities of the perturbation analysis are designed to be constant relative to $ϵ$. Or if you do it, then $δ$ is a constant and $1+δϵ$ is a linear, not a quadratic perturbation of $1$. By the claim of the task, you will not find the region of divergent resonance with this linear approach.
The perturbation system for the quadratic approach is
\begin{align}
\partial_{\tau\tau}x_0+x_0&=0, \tag{1}
\\
\partial_{\tau\tau}x_1+x_1&=-x_0\cos t, \tag{2}
\\
\partial_{\tau\tau}x_2+x_2&=-2\partial_{T\tau}x_0-cx_0-x_1 \cos t, \tag3
\end{align}
New version using standard trigonometric parametrization
The solution of the base equation has, as you also used, the form
\begin{align}
x_0(τ,T) &= A(T)\cos(τ)+B(T)\sin(τ)
\end{align}
Then the next order term satisfies
\begin{align}
∂_{ττ}x_1+x_1&=-\frac{A(T)}2(1+\cos2τ)-\frac{B(T)}2\sin2τ\\
x_1(τ,T) &= -\frac{A(T)}2 + \frac{A(T)}6\cos2τ + \frac{B(T)}6\sin2τ
\end{align}
Note that the frequency one components of the homogeneous/complementary solution were left out, as they would only replicate some fraction of the base solution.
In the next perturbation term we encounter the first $T$ derivative and the frequency perturbation, as well as the twice shifted base frequency, all sources for resonance.
\begin{align}
∂_{ττ}x_2+x_2&=-2∂_{Tτ}x_0 - cx_0 - \cosτ\,x_1 \\
&= 2A'(T)\sinτ - 2B'(T)\cosτ -cA(T)\cosτ - cB(T)\sinτ\\&\qquad
+ \frac{A(T)}2\cosτ - \frac{A(T)}{12} (\cos3τ+\cosτ) - \frac{B(T)}{12}(\sin3τ + \sinτ)
\end{align}
To avoid resonance, the frequency one terms have to cancel,
\begin{align}
2A'(T)&=cB(T) + \frac{B(T)}{12}&&=\left(\frac1{12}+c\right)B(T)
\\
2B'(T)&=-cA(T) + \frac{A(T)}2 - \frac{A(T)}{12}&&=\left(\frac{5}{12}-c\right)A(T)
\end{align}
This has a pair of real eigenvalues if $\left(\frac1{12}+c\right)\left(\frac{5}{12}-c\right)>0$, that is, for $-\frac1{12}<c<\frac{5}{12}$. As one of them is positive, this gives an exponentially growing term in the solution, leading to divergence as per the claim.
Old version using phase parametrization
An equivalent form of the solution of the base equation is
\begin{align}
x_0(τ,T)&=A(T)\cos(τ+B(T)), \\
x_0(τ,T)\cos(τ)&=\frac{A(T)}2[\cos(2τ+B(T))+\cos(B(T))],
\end{align}
Then the first perturbation term is computed as
\begin{align}
x_1(τ,T)&=\frac{A(T)}6\cos(2τ+B(T))-\frac{A(T)}2\cos(B(T)),
\\
x_1(τ,T)\cos(τ)&=\frac{A(T)}{12}[\cos(3τ+B(T))+\cos(τ+B(T))]\\
&~~~~-\frac{A(T)}2[\cos(τ+B(T))\cos^2(B(T))+\sin(τ+B(T))\sin(B(T))\cos(B(T))]
\\[.5em]
&=\frac{A(T)}{12}\cos(3τ+B(T))-\frac{A(T)}{12}[2+3\cos(2B(T))]\cos(τ+B(T))
\\&~~~~-\frac{A(T)}4\sin(2B(T))\sin(τ+B(T)),
\end{align}
Note that the frequency one components of the homogeneous/complementary solution were left out, as they would only replicate some variant of the base solution.
In the equation for the next term we encounter resonance which has to be resolved using the dependence of the parameters $A(T)$ and $B(T)$ on $T$.
\begin{align}
∂_{ττ}x_2+x_2&=2\Bigl[A'(T)\sin(τ+B(T))+A(T)B'(T)\cos(τ+B(T))\Bigr]\\
&~~~~-cA(T)\cos(τ+B(T))-x_1(τ,T)\cos(τ)
\\[.5em]
&=\Bigl[2A'(T)+\frac{A(T)}4\sin(2B(T))\Bigr]\sin(τ+B(T))\\
&~~~~+\Bigl[2A(T)B'(T)-cA(T)+\frac{A(T)}{12}[2+3\cos(2B(T))]\Bigr]\cos(τ+B(T))\\
&~~~~-\frac{A(T)}{12}\cos(3τ+B(T))
\end{align}
To remove the explicit growth in the form of polynomial coefficients we need to set the coefficients of the first two resonance terms, that is, of $\cos(τ+B(T))$ and $\sin(τ+B(T))$, to zero, giving the system
\begin{align}
8A'(T)&=-A(T)\sin(2B(T))\\
8B'(T)&=4c-\frac23-\cos(2B(T))
\\[1em]
\implies A(T)&=\frac{C}{4c-\frac23-\cos(2B(T))}
\end{align}
Now as long as $4c-\frac23\in(-1,1)\iff c\in(-\frac1{12},\frac{5}{12})$, the second equation has an infinite amount of simple equilibria, keeping $B(T)$ bounded between two of them. Asymptotically $B(T)$ will converge to the next stable equilibrium, that is, where $\sin(2B(T))$ is negative. But then the first equation will have a positive growth factor, so that the amplitude $A(T)$ is exponentially growing.
Outside that range, $A(T)$ is a periodic function of $B(T)$, and thus bounded.
This second approach is visibly more complicated due to multiple different applications of trigonometric identities, than the first one, and much harder to check for errors.
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$.
Best Answer
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.