# solving coupled differential equations

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 :

• Your initial conditions should be enclosed within your curly braces with the equations. I reached the maximum number of steps at t=0.096, so you will have to investigate that. Commented Jan 9, 2020 at 11:40
• Thank you. It worked. Commented Jan 10, 2020 at 5:45

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]}


• Thank you. What I need to do to plot the (2,2) element of the matrix M? Commented Jan 10, 2020 at 5:44
• Use 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}, Flatten[{Variables[M], Variables[M2]}], {t, 0, 10^5}, MaxSteps -> Infinity] and then {Plot[Variables[M][[4]] /. First[s] // Re, {t, 0, 10^5}], Plot[Variables[M][[4]] /. First[s] // Im, {t, 0, 10^5}]} Commented Jan 10, 2020 at 11:24