Solve Mathieu Equation $y''+(1+\epsilon \cdot \sin(3t))y=0, \;y(0)=u, \;y'(0)=v$ I just need to know the solution to this differential equation. I tried wolfram alpha and some other online resources but none worked. 
$$y''(t)+(1+\epsilon \cdot \sin(3t))y(t)=0, \; y(0)=u, \; y'(0)=v$$
 A: $$\frac{d^2y}{dt^2}+(1+\epsilon \sin(3t))y(t)=0 \tag 1$$
$t=\frac{2}{3}\tau \quad\to\quad \frac{d^2y}{dt^2}=\frac{9}{4}\frac{d^2y}{d\tau^2}$ 
$$\frac{d^2y}{d\tau^2}+\left(\frac{4}{9}+\frac{4}{9}\epsilon \sin(2\tau)\right)y(\tau)=0$$
$\tau=x+\frac{\pi}{4} \quad\to\quad \frac{d^2y}{d\tau^2}=\frac{d^2y}{dx^2}  \quad $and$\quad \sin(2\tau)=\cos(2x)$
$$\frac{d^2y}{dx^2}+\left(\frac{4}{9}+\frac{4}{9}\epsilon \cos(2x)\right)y(x)=0$$
This is the Mathieu ODE : $\quad\frac{d^2y}{dx^2}+\left(a-2b \cos(2x)\right)y(x)=0\quad$ with $\begin{cases}a=\frac{4}{9}\\b=-\frac{2\epsilon}{9} \end{cases}$
The general solution involves the Mathieu functions :
$y(x)=c_1\text{MathieuC}(a,b,x)+c_2\text{MathieuS}(a,b,x)$
$$y(x)=c_1\text{MathieuC}\left(\frac{4}{9},-\frac{2\epsilon}{9},x\right)+
c_2\text{MathieuS}\left(\frac{4}{9},-\frac{2\epsilon}{9},x\right)$$
$$y(\tau)=c_1\text{MathieuC}\left(\frac{4}{9},-\frac{2\epsilon}{9},\tau-\frac{\pi}{4}\right)+c_2\text{MathieuS}\left(\frac{4}{9},-\frac{2\epsilon}{9},\tau-\frac{\pi}{4}\right)$$
The general solution of the ODE $(1)$ is :
$$y(t)=c_1\text{MathieuC}\left(\frac{4}{9}\:,\:-\frac{2\epsilon}{9}\:,\:\frac{3}{2}t -\frac{\pi}{4}\right)+c_2\text{MathieuS}\left(\frac{4}{9}\:,\:-\frac{2\epsilon}{9}\:,\:\frac{3}{2}t -\frac{\pi}{4}\right)$$
The coefficients $c_1$ and $c_2$ have to be determined according to the conditions $y(0)=u$ and $y'(0)=v$. This is theoeticaly possible, but in fact very arduous, due to the antiderivatives of the Mathieu functions which have no simple closed form.
NOTE :
If the real question is not to solve Eq.$(1)$ , but to find an approximate solution for small $\epsilon$ , the problem is very different. Using the exact solution in terms of Mathieu function is not a smart method. It would be a bad way, very complicated. 
It is much simpler to expend $y(\epsilon,t)$ in power series of $\epsilon$ :
$$y(t)=f_0(t)+\epsilon f_1(t)+\epsilon^2 f_2(t)+...$$
Putting it into Eq.$(1)$ :
$(f_0''+\epsilon f_1''+\epsilon^2 f_2''+...)+(1+\epsilon \sin(3t))(f_0+\epsilon f_1+\epsilon^2 f_2+...)=0$
$(f_0''+f_0)+\epsilon \left(f_1''+f_1+\sin(3t)f_0\right)+\epsilon^2 \left(f_2''+f_2+\sin(3t)f_1\right)+...=0$
The first approximate is $y(t)\simeq f_0(t)$ with $\epsilon=0$ :
$f_0''(t)+f_0(t)=0$ with conditions $f_0(0)=u$ and $f_0'(0)=v\quad\to\quad $
$$f_0(t)=u\:\cos(t)+v\:\sin(t)$$
For the second approximate : $\quad f_1''+f_1+\sin(3t)f_0=0\quad$ with conditions $f_1(0)=f'_1(0)=0$ since the boundary conditions are already satisfied above.
$f_1''+f_1+\sin(3t)(u\:\cos(t)+v\:\sin(t))=0$
$$f_1(t)=-\frac{4}{15}\sin^3(t)\left(11u\cos{\frac{t}{2}}+3u\cos{\frac{3t}{2}}+u\cos{\frac{5t}{2}}+ v\sin{\frac{t}{2}}+3v\sin{\frac{3t}{2}}+v\sin{\frac{5t}{2}}\right)$$
$$y(t)\simeq u\:\cos(t)+v\:\sin(t) +\epsilon \:f_1(t)$$
For even better approximate, one could continue on the same way with $f_2''+f_2+\sin(3t)f_1=0$ and conditions $f_2(0)=f_2'(0)=0$.
