show/hide this revision's text 3 added 10 characters in body

Another thing to consider: introduce new independent variable $\tau=sin(\alpha \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.
show/hide this revision's text 2 added 98 characters in body; deleted 46 characters in body

Another thing to consider: introduce new independent variable $\tau=sin(\alpha z+\Phi_a)$. Also note that we have $$cos(\alpha z+\Phi_b)=cos(\alpha z+\Phi_a+(\Phi_b-\Phi_a))$$ z+\Phi_b)=cos(x+\delta)$$ and $$ cos(\alpha z+\Phi_a+(\Phi_b-\Phi_a))=cos(\alpha z+\Phi_a)cos(\Phi_b-\Phi_a)-sin(\alpha z+\Phi_a)sin(\Phi_b-\Phi_a).$$ 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(\Phi_b-\Phi_a)-\alpha\tau sin(\Phi_b-\Phi_a)/\sqrt{1-\alpha^2\tau^2}\rightleft(cos(\delta)-\displaystyle\frac{\alpha\tau}{\sqrt{1-\alpha^2\tau^2}} sin(\delta)\right) y_2, y_2,$$ while the second one is simply $$ dy_2/d\tau=A y_1$$ which y_1. $$ This new system might be somewhat easier to investigate, be it analyticlaly analytically or numerically.
show/hide this revision's text 1

Another thing to consider: introduce new independent variable $\tau=sin(\alpha z+\Phi_a)$. Also note that we have $$cos(\alpha z+\Phi_b)=cos(\alpha z+\Phi_a+(\Phi_b-\Phi_a))$$ and $$ cos(\alpha z+\Phi_a+(\Phi_b-\Phi_a))=cos(\alpha z+\Phi_a)cos(\Phi_b-\Phi_a)-sin(\alpha z+\Phi_a)sin(\Phi_b-\Phi_a).$$ Then your system becomes $$ dy_1/d\tau=B \left(cos(\Phi_b-\Phi_a)-\alpha\tau sin(\Phi_b-\Phi_a)/\sqrt{1-\alpha^2\tau^2}\right) y_2, dy_2/d\tau=A y_1$$ which might be somewhat easier to investigate, be it analyticlaly or numerically.