Method of multiple scales on Mathieu's equation I encountered a problem in Strogatz's "Nonlinear Dynamics and Chaos." Specifically 7.6.18. He takes the following equation : $$\ddot{x}+(a+\epsilon \cos t)x=0$$ where a is close to 1. There he asks us to use two timing method with slow time $\epsilon^2t$
and show the following:

For $1-\frac{\epsilon^2}{12}+o(\epsilon^4)<a<1+\frac{5\epsilon^2}{12}+o(\epsilon^4)$, the solution is unbounded.

I set $a=1+\epsilon\delta$ where $\delta$ is a $O(1)$ constant so the equation turns to $\ddot{x}+x+\epsilon(\delta+\cos t)x=0$.
I take $x(t,\epsilon)=x_0(\tau,T)+\epsilon x_1(\tau,T)+\epsilon^2x_2(\tau,T)+O(\epsilon^3)$, where $\tau=t$ and $T=\epsilon^2 t$
I get 
\begin{align}
\frac{dx}{dt}&=\frac{\partial x}{\partial \tau}+\epsilon^2\frac{\partial x}{\partial T}
\\
\frac{d^2x}{dt^2}
&=\partial_{\tau\tau}x+2\epsilon^2\partial_{T\tau }x+O(\epsilon^4)\\
&=\partial_{\tau\tau}(x_0+\epsilon x_1+\epsilon^2x_2+O(\epsilon^3))+2\epsilon^2\partial_{T\tau}(x_0+\epsilon x_1+..)+O(\epsilon^3)\\
&=\partial_{\tau\tau}x_0+\epsilon\partial_{\tau\tau}x_1+\epsilon^2(\partial_{\tau\tau}x_2+2\partial_{T\tau}x_0)+O(\epsilon^3)
\end{align}
so $\ddot x+x+\epsilon(\delta+\cos t)x=0$ gets transformed into:
$$
\partial_{\tau\tau}x_0+x_0+\epsilon(\partial_{\tau\tau}x_1+x_1+x_0(\delta+\cos t))+\epsilon^2(\partial_{\tau\tau}x_2+2\partial_{T\tau}x_0+x_1(\delta+\cos t))=O(\epsilon^3)
$$
So the perturbation equations till 2nd order are:
\begin{align}
\partial_{\tau\tau}x_0+x_0&=0 \tag{1}
\\
\partial_{\tau\tau}x_1+x_1&=-x_0(\delta+\cos t)\tag{2}
\\
\partial_{\tau\tau}x_2+x_2&=-2\partial_{T\tau}x_0-x_1(\delta+\cos t)\tag3
\end{align}
general solution to equation 1 is 
$$
x_0(\tau,T)=A(T)\cos\tau+B(T)\sin\tau\tag4
$$
using 4 in 2 we get:
$$
\partial_{\tau\tau}x_1+x_1=-A\delta \cos\tau-B\delta \sin\tau-\frac{A}{2}(\cos2\tau+1)-\frac{B}{2}\sin2\tau
$$
for non secular solutions in equation 2, we set coefficients of resonant terms to 0 giving $A=0$ and $B=0$ which in turn gives us $x_0=0$
and equation 2 becomes $\partial_{\tau\tau}x_1+x_1=0$ and equation 3 becomes $\partial_{\tau\tau}x_2+x_2=-x_1(\delta+\cos t)$ 
This is what we started with for $x_0$ so $x_1=0$ and so on.
So just doing the two timing with $T=\epsilon^2t$ does not seem to work.
I tried using three timescales with $x(t,\epsilon)=x_0(\tau,T,\sigma)+\epsilon x_1(\tau,T,\sigma)+....$ where $\tau=t$, $T=\epsilon t$, $\sigma=\epsilon^2 t$
but the algebra got a bit unwieldy and suggested that this was not what the author wanted us to do.
Can anybody prove the above result?
Any suggestions are appreciated. 
 A: 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.
