Homogeneous linear differential equation system with simple periodical coefficient matrix Hello, I encountered the following system of linear first-order differential equations:
$y'(z)=A(z) y(z)$
where
$y(z): R \rightarrow R^2$ and
$A(z)=\begin{pmatrix} 
0 & B Cos(\alpha z + \Phi_b) \cr 
A Cos(\alpha z + \Phi_a) & 0 \end{pmatrix}$
All quantities are real and constant if not explicitly variable. (Probably - if anyone wants to talk about approximations, which would be ok to some extent -) $A$,$B$ are $<1$, while $\alpha>10^3$ and $z<10^{-2}$.
I found Floquet's theorem, which only gives my the structure of the solution but isn't of much help for me. 
In approximation, I thought about using the Magnus expansion.
I added math.SP, as I got the suggestion to transform the equation to a Schrodinger equation with $Y(z)=A_{1,2}(z)^{1/2} y(z)$ but that did not produce any new ideas for a solution.
I am at the end my (physicist's) knowledge, hope anyone has an idea.
 A: OK. First of all, change $(y_1,y_2)$ to $(y_1,\sqrt{B/A}y_2)$ and the time $z$ to $t=\alpha z+\frac{\Phi_A+\Phi_B}{2}$. Then we'll get the system with the matrix 
$$
A(t)=\begin{pmatrix} 
0 & C\cos (t + \Phi) \cr 
C \cos(t- \Phi) & 0 \end{pmatrix}
$$
where $C=\alpha^{-1}\sqrt{AB}$ and $\Phi=(\Phi_B-\Phi_A)/2$. We want our approximation to be decent on $[\alpha z_{\min} ,\alpha z_{\max}]\subset [-\alpha,\alpha]$. Noting that $\cos(t-\Phi)=\cos(\Phi-t)$, we see that $A(t)=A^*(2\pi-t)$, which immediately tells us that the monodromy matrix $M$ from $0$ to $2\pi$ is self-adjoint. Also, denoting $\psi(t)=C\cos(\Phi+t)$=ce^{it}+\bar c{e^{-it}}$ with $c=\frac 12Ce^{i\Phi}, we see that the fundamental matrix $M(t)$ of the solution on $[0,2\pi]$ can be obtained (by the standard Piquard iterations) as the sum
$$
M(t)=\begin{pmatrix} 
1 & 0 \cr 
0 & 1 \end{pmatrix}+
\begin{pmatrix} 
0 & \psi_1(t) \cr 
-\psi_1(-t) & 0 \end{pmatrix}
$$
$$+
\begin{pmatrix} 
 -\psi_2(t) & 0 \cr 
0 & -\psi_2(-t)  \end{pmatrix}+
\begin{pmatrix} 
0 & -\psi_3(t) \cr 
\psi_3(-t) & 0 \end{pmatrix}+
O(C^4)
$$ 
Where $\psi_0(t)=1$ and $\psi_{k+1}(t)=\int_0^t\psi(s)\psi_k(-s)ds$.
You can write a long series but I want to convince you that 4 first terms are enough for your problem. We can find the first two $\psi$'s:
$$
\psi_1=\frac 1i[(ce^{it}-\bar c{e^{-it}})-(c-\bar c)
$$
$$
\psi_2=-\frac{c-\bar c}{i}\psi_1+\frac 1i(c^2-\bar c^2)t+\frac{|c|^2}2(e^{2it}+e^{-2it})-|c|^2
$$
and the linear term in $\psi_3$, which is
$$
(c^2-\bar c^2)t[(c-\bar c)+ce^{it}-\bar ce^{-it}]
$$
Plugging in $t=2\pi$, we see that the monodromy matrix is $\begin{pmatrix} 
1-2\pi v & 4\pi vs \cr 
4\pi vs & 1+2\pi v \end{pmatrix}+\begin{pmatrix} 
O(C^4) & O(C^5) \cr 
O(C^5) & O(C^4)\end{pmatrix}$
with $v=2\Im (c^2)=\frac 12 C^2\sin 2\Phi= \frac 12 \alpha^{-2}AB\sin(\Phi_B-\Phi_A)$, $s=c-\bar c=C\sin\Phi=\alpha^{-1}\sqrt{AB}\sin\frac{\Phi_B-\Phi_A}2$. 
Now, the life is easy: the growth/decay part is essentially given by the matrix 
$$
G(t)=\begin{pmatrix} 
e^{-vt} & 0 \cr 
0 & e^{vt} \end{pmatrix}
$$ 
the rotation part is essentially given by 
$$
T=\begin{pmatrix} 
1 & -s \cr 
s & 1 \end{pmatrix}
$$
and the oscillation is essentially given by 
$$
H(t)=\begin{pmatrix} 
1-\widetilde\psi_2(t) & \psi_1(t) \cr 
-\psi_1(-t) & 1-\widetilde\psi_2(-t) \end{pmatrix}
$$
where $\widetilde \psi_2$ is $\psi_2$ with the term $\frac 1i(c^2-\bar c^2)t$ removed. 
So, $M(t)\approx H(t)T^{-1}G(t)T$ (well, $T^{-1}GT$ and $H$ don't commute but the commutator effect is of size $C^3$).
This should work just fine letting you to see just enough in your range. 
P.S. My original answer had an error in that I neglected the rotation of eigenvectors, but now it should be fine even for fairly large $C$ like $0.1$. Check agaist your numerics and see if it works well enough for you.
A: Set $\alpha z=w$, you get the new system
$$dy/dw={1\over \alpha}\pmatrix{0&B\cos(w+\Phi_b)\cr A\cos(w+\Phi_a)&0}y.$$
Since you are interested in a case where $\alpha$ is large and $w$ is of moderate size, you can try an expansion of the solution in powers of $1/\alpha$.
A: Another thing to consider: introduce new independent variable $\tau=(1/\alpha)sin(\alpha z+\Phi_a)$. Also note that we have
$$cos(\alpha z+\Phi_b)=cos(x+\delta)$$
and
$$
cos(x+\delta)=cos(x)cos(\delta)-sin(x)sin(\delta).$$
where $x=\alpha z+\Phi_a$ and $\delta=\Phi_b-\Phi_a$.
Then the first equation of your system becomes
$$
dy_1/d\tau=B \left(cos(\delta)-\displaystyle\frac{\alpha\tau}{\sqrt{1-\alpha^2\tau^2}} sin(\delta)\right) y_2,$$ 
while the second one is simply 
$$
dy_2/d\tau=A y_1.
$$
This new system might be somewhat easier to investigate, be it analytically or numerically.
