The Mathieu equation is a second-order linear differential equation given by $$y''(t) + [a - 2q\cos(2t)]y(t) = 0$$ There are two special functions defined as linearly independent solutions to Mathieu's equation:
1. The even Mathieu cosine function which I will denote as $MC(a,q,t)$, satisfying the initial conditions $MC(a,q,0)=1$ and $MC'(a,q,0)=0$ (where prime denotes $t$ derivative).
2. The odd Mathieu sine function, $MS(a,q,t)$, satisfying $MS(a,q,0)=0$ and $MS'(a,q,0)=1$.
For $q=0$, Mathieu's equation reduces down to the familiar differential equation $$y''(t) + ay(t) = 0 \tag{*}$$ for which we have even and odd solutions $\cos(\sqrt{a}t)$ and $\sin(\sqrt{a}t)$. Let $\epsilon = \frac{2q}{a}$ and suppose that $2q \ll a$ so that $\epsilon \ll 1$. Therefore Mathieu's equation can be written as $$\frac{1}{a}y''(t) + [1-\epsilon \cos(2t)]y(t) = 0$$ where we view it as a perturbed form of equation $(*)$. It follows that for the case of small $\epsilon$, the functions $MS$ and $MC$ should be approximated by $\sin$ and $\cos$ plus small perturbative terms.
The problem with such a perturbative solution is the appearance of secular terms at $\mathcal{O}(\epsilon^2)$, i.e. terms of the form $t\cos(\sqrt{a} t)$ or $t\sin(\sqrt{a} t)$ which diverge with large $t$, whereas the functions $MC$ and $MS$ appear to be uniformly bounded.
I have tried several standard methods for obtaining uniform perturbative solutions, i.e. the Poincaré - Lindstedt method and multiple scale analysis but neither method successfully gets rid of the secular terms. Does anyone know how I can (or if it's even possible) to obtain an approximation for $MC$ and $MS$ which is reasonably accurate for all $t$?
It is possible to get a leading-order approximation to the solution via multiple scales relatively easily, as follows. First, for convenience, define $\tau = a^{1/2}t$, such that the equation reads $$y''(\tau)+(1-\epsilon \cos(2a^{-1/2}\tau))y(\tau)=0. $$
Since the unperturbed problem has a constant frequency, we might suspect a multiple scales approach to work. Define a slow time scale $T=\epsilon \tau$ and consider $y\equiv y(\tau,T)$, such that $\frac{\mathrm{d}}{\mathrm{d}\tau} = \partial_\tau + \epsilon \partial_T$. The ODE becomes $$ (\partial_\tau+\epsilon\partial_T)^2 y(\tau,T) + (1-\epsilon\cos(2a^{-1/2}\tau))y(\tau,T)=0. $$ Expand the unknown function $y(\tau,T)=y_o(\tau,T)+\epsilon y_1(\tau,T)+o(\epsilon)$. At order one, we find the unperturbed problem $$\partial_\tau^2y_o+y_o=0, $$ whose solution is $y_o(\tau,T)=A(T)e^{i\tau}+c.c.$, where $c.c.$ means "complex conjugate". At order $\epsilon$, $$ \partial_\tau^2y_1+y_1 = -2 \partial_\tau \partial_T y_o + \cos(2a ^{-1/2}\tau)y_o. $$ The RHS is $$-2\partial_\tau \partial_T y_o + \cos(2a ^{-1/2}\tau)y_o = -2(iA_T e^{i\tau} +c.c.)+\frac{1}{2} \left(e^{2a^{-1/2}i\tau}+e^{-2a^{-1/2}i\tau} \right) \left( A e^{i\tau} +c.c.\right) $$ There are two cases to be distinguished at this point. If $|2a^{-1/2}\pm 1| \neq 1$, then there are no resonances stemming from the product term. In this case, the solvability condition gives simply $$A_T =0, $$ s.t. the leading-order result is $y_o(\tau,T)\equiv y_o(\tau)= A e^{i\tau} +c.c.=A e^{i\sqrt{a}t} +c.c.$, with $A\in \mathbb{C}$ a constant. If $a=1$, however, then there's a resonant contribution from the product term and the solvability condition reads $$ \partial_T A = - \frac{1}{4}i A^* $$ which, upon differentiating again w.r.t. $T$, yields $\partial_T^2 A = \frac{1}{16} A $, which is easily solved by $A(T) =(\alpha-i\beta)e^{T/4}+(\gamma-i\delta)e^{-T/4}$ for real constants $\alpha,\beta,\gamma,\delta$. Plugging this into the zeroth order solution gives $$y_o(\tau,T)= \left(\alpha e^{T/4}+\gamma e^{-T/4}\right) \cos(\tau) + \left(\beta e^{T/4} + \delta e^{-T/4}\right)\sin(\tau).$$ Using that in this case $t=\tau$, by definition $T=\epsilon \tau$ back in, $$ y_o(t,T)= \left(\alpha e^{\epsilon t/4}+\gamma e^{-\epsilon t/4}\right) \cos(t) + \left(\beta e^{\epsilon t/4} + \delta e^{-\epsilon t/4}\right)\sin(t). $$ Both approximations at $a=1$ or $a\neq 1$ are valid until $t=O(1/\epsilon)$. If higher order accuracy is needed, one can introduce a third time-scale $\tilde{T}=\epsilon^2 t$ and go one order further.