solving coupled differential equations

Question

I am solving 6 differential equations in matrix form as:

\[Theta] = 34.5*\[Pi]/180;
m2 = 1/2*7.3*10^-5*{{Cos[
     2 \[Theta]], -Sin[2 \[Theta]]}, {-Sin[2 \[Theta]], -Cos[
      2 \[Theta]]}}
M = Array[Subscript[Subscript[\[Rho], 1], #1, #2][t] &, {2, 2}]
M2 = Array[Subscript[c, #1, #2][t] &, {2, 2}]
M0 = {{1, 0.00001}, {0.00001, 0}};
ci = Thread[Flatten[M] == Flatten[M0]] /. {t -> 0};
M20 = {{2, 0.00001}, {0.00001, 0}};
di = Thread[Flatten[M2] == Flatten[M20]] /. {t -> 0};

s = NDSolve[{ 


 I D[M, t] == (m2/10 + 10^5*2*M2).M - M.(m2/10 + 10^5*2*M2), 
   I D[M2, t] == (m2/11 + 10^5*2*M).M2 - M2.(m2/11 + 10^5*2*M)}, ci, 
  di, {Variables[M], Variables[M2]}, {t, 0, 10^7}]

But it is showing the error as :

Subscript[c, 1,1][0]==2 cannot be used as a variable. >> Can someone please help me with this?

Answer 1

We can normalize the equations by 10^5. So that the oscillation period is of the order of 1, we make the change t1->t*10^5. Then it is possible to integrate up to 10^5, since the dynamics are also visible on this scale.

\[Theta] = 34.5*\[Pi]/180;
m2 = 1/2*7.3*10^-5*{{Cos[
      2 \[Theta]], -Sin[2 \[Theta]]}, {-Sin[2 \[Theta]], -Cos[
       2 \[Theta]]}};
M = Array[Subscript[Subscript[\[Rho], 1], #1, #2][t] &, {2, 2}];
M2 = Array[Subscript[c, #1, #2][t] &, {2, 2}];
M0 = {{1, 0.00001}, {0.00001, 0}};
ci = Thread[Flatten[M] == Flatten[M0]] /. {t -> 0};
M20 = {{2, 0.00001}, {0.00001, 0}};
di = Thread[Flatten[M2] == Flatten[M20]] /. {t -> 0};

s = NDSolve[{I D[M, t] == (m2/10^6 + 2*M2).M - M.(m2/10^6 + 2*M2), 
   I D[M2, t] == (m2/11/10^5 + 2*M).M2 - M2.(m2/11/10^5 + 2*M), ci, 
   di}, {Variables[M], Variables[M2]}, {t, 0, 10^5}, 
  MaxSteps -> Infinity]

Two-scale visualization

    {Plot[Evaluate[Variables[M] /. First[s][[1]] // Im], {t, 0, 10^1}, 
  PlotRange -> All, Frame -> True, Axes -> False, 
  PlotLegends -> Variables[M]], 
 Plot[Evaluate[Variables[M2] /. First[s][[2]] // Im], {t, 0, 10^1}, 
  PlotRange -> All, Frame -> True, Axes -> False, 
  PlotLegends -> Variables[M2]]}

{Plot[Variables[M] /. First[s][[1]] // Im, {t, 0, 10^5}, 
  PlotRange -> All, Frame -> True, Axes -> False], 
 Plot[Variables[M2] /. First[s][[2]] // Im, {t, 0, 10^5}, 
  PlotRange -> All, Frame -> True, Axes -> False]}