I'm trying to solve the driven Mathieu equation $x''+\beta x'+(a-2q\cos{\Omega t})\frac{\Omega^2}{4}x=f(t)$ for both zero and non-zero $\beta$.

I can write down an analytic solution using the variation of parameters method (writing $x(t)=u_1(t)x_1(t)+u_2(t)x_2(t)$ where $x_1(t)$ and $x_2(t)$ are independent solutions of the un-driven Mathieu equation, leading to $u_1(t)=-\int^t_\alpha\frac{x_2(\tau)f(\tau)}{W(\tau)}{\rm d}\tau$ and similar for $u_2(t)$), however that produces intractable integrations. What's the best method for (approximately) solving this ODE?

So far, I'm using approximate expressions for my solutions, $x_1(t)=\cos\omega t(1-\frac{q}{2}\cos\Omega t)$ and $x_2(t)=\sin\omega t(1-\frac{q}{2}\cos\Omega t)$ ($\omega=\frac{\Omega}{2}\sqrt{a+\frac{q^2}{2}}$,$q<<1$, $\beta=0$), and then throwing terms away by order of highest powers of $q$ in my integration until Mathematica can handle it. For the simple case $f(t)=\rm const$ I have to go all the way to $W(t)=\omega$ before the integral is doable.

Is this in principle the best method of solution or am I missing a trick (I'm an experimental physicist, so my mathematical skills can be somewhat lacking)? How can I assess how accurate my answer is, beyond comparing it visually with a numerical integration?