function dydt=equations(t,y,ky,kz,B0z)

%y(1)=D, y(2)=dD/dt, y(3)=bx, y(4)=dbx/dt, y(5)=by, y(6)=dby/dt

q=1.5;
kx=q*ky*t;

dydt=[ y(2)
     
      -kz*kz*y(1)+kx*kz*y(3)+ky*kz*y(5)
      
       y(4)
       
      kx*kz*y(1)-(kx^2+(kx*B0z)^2+(kz*B0z)^2-2.0*q)*y(3)-kx*ky*(1+B0z^2)*y(5)+2.0*y(6)
      
       y(6)
       
      ky*kz*y(1)-kx*ky*(1+B0z^2)*y(3)-2.0*y(4)-(ky^2+(ky*B0z)^2+(kz*B0z)^2)*y(5)];
       
 




