s = SessionTime[];

\(\(	
x0 = .2; y0 = .3; k = 5.; \)\)

\(\(
\(While[
      True, {blst = {}; lst = {}; 

        Do[{
			lst = {}; 
				
            Do[
					{
                mathieu = {\(x'\)[t] == y[t], 
                    \(y'\)[t] == \(-\((s^2 + e\ Cos[2  π\ t])\)\) x[t], 
                    x[0] == .2, y[0] == .3}; 	

                sol = NDSolve[mathieu, {x[t], y[t]}, {t, 0, 2}]; 

                x1 = sol[\([1, 1, 2]\)] /. {t -> 1}; 

                x2 = sol[\([1, 1, 2]\)] /. {t -> 2}; 

                y1 = sol[\([1, 2, 2]\)] /. {t -> 1}; 

                y2 = sol[\([1, 2, 2]\)] /. {t -> 2}; 

                m = Solve[{a\ x0 + b\ y0 == x1, c\ x0 + d\ y0 == y1, 
                      a\ x1 + b\ y1 == x2, c\ x1 + d\ y1 == y2}]; 

                A = {{m[\([1, 1, 2]\)], m[\([1, 2, 2]\)]}, {m[\([1, 3, 2]\)], 
                      m[\([1, 4, 2]\)]}}; 
eig = Eigenvalues[A]; 

                						
                lst = Append[lst, 
                    If[Chop[Abs[eig[\([1]\)]] - Abs[eig[\([2]\)]]] == 0, 
                      Blue, Red]]}, {e, \(-50\), 50, 10  k}]; 
		
            blst = Append[blst, lst]}, {s, 0, 10, k}]; 

        Show[Graphics[RasterArray[Transpose[blst]], 
            PlotLabel -> Ceiling[\((\((10/k)\) + 1)\)]]]; 

        Print[Ceiling[\((SessionTime[] - s)\)], " sec."]; 
	
        s = SessionTime[]; k = k/4}]; \)
\)\)