Matrix ODE gives badly-behaved numerical solution [duplicate]

Question

I am trying to solve matrix ODE:

(*change global precision*)
$PreRead = (# /.s_String /; StringMatchQ[s, NumberString] && Precision@ToExpression@s == MachinePrecision :> s <> "`44." &);
(*define matrix components*)
ρ[t_] := Table[a[i, j][t], {i, 4}, {j, 4}] 
slnn = Flatten[Table[a[i, j][t], {i, 4}, {j, 4}]];

(*define initial condition for ODE*)
init = ( {{0, 0, 0, 0},{0, 1, 0, 0},{0, 0, 0, 0},{0, 0, 0, 0}} ); γ = 1; θ = Pi/2;
(*now define constants*)
r[i_, j_] := (i - j) di; di = 1/2*10^-7; μ = 10^-10; ω0 =2*10^15; c = 3*10^8; N[2 Pi/(ω0/c)]
σ1p = ( {{0, 0, 1, 0},{0, 0, 0, 1},{0, 0, 0, 0},{0, 0, 0, 0}} ); σ1m = ( {{0, 0, 0, 0},{0, 0, 0, 0},{1, 0, 0, 0},{0, 1, 0, 0}} ); σ2p = ( {{0, 1, 0, 0},{0, 0, 0, 0},{0, 0, 0, 1},{0, 0, 0, 0}} ); σ2m = ( {{0, 0, 0, 0},{1, 0, 0, 0},{0, 0, 0, 0},{0, 0, 1, 0}} );
Λ12 = 3/4 (24 Cos[1/3] + 9 Sin[1/3]);
Λ21 = 3/4 (24 Cos[1/3] + 9 Sin[1/3]);
sh = Sinh[2]; ch = Cosh[2];
Λp21 = -0.47260812;
Λp11 = -0.47260812;
(*now we solve ODE*)
solsq = NDSolve[{ρ'[
      t] == -sh ch (((1/2 + 
              I Λp21) (σ2p.σ1p.ρ[
                t] + ρ[t].σ2p.σ1p) + (1/2 + 
              I Λp21) (σ1p.σ2p.ρ[
                t] + ρ[t].σ1p.σ2p)) + ((1/2 - 
              I Λp21) (σ2m.σ1m.ρ[
                t] + ρ[t].σ2m.σ1m) + (1/2 - 
              I Λp21) (σ1m.σ2m.ρ[
                t] + ρ[t].σ1m.σ2m)) - ((1 + 
              2 I Λp21) σ2p.ρ[
              t].σ1p + (1 + 
              2 I Λp21) σ1p.ρ[
              t].σ2p + (1 + 
              2 I Λp11) σ1p.ρ[
              t].σ1p + (1 + 
              2 I Λp11) σ2p.ρ[
              t].σ2p) - ((1 - 
              2 I Λp21) σ1m.ρ[
              t].σ2m + (1 - 
              2 I Λp21) σ2m.ρ[
              t].σ1m + (1 - 
              2 I Λp11) σ1m.ρ[
              t].σ1m + (1 - 
              2 I Λp11) σ2m.ρ[
              t].σ2m)) - 
      1/2 (sh + 
         1) ((ρ[t].σ1p.σ1m + ρ[
             t].σ1p.σ2m + ρ[
             t].σ2p.σ1m + ρ[
             t].σ2p.σ2m) + \
(σ1p.σ1m.ρ[t] + σ1p.σ2m.ρ[
             t] + σ2p.σ1m.ρ[
             t] + σ2p.σ2m.ρ[t]) - 
     2 (σ1m.ρ[t].σ1p + σ1m.ρ[
          t].σ2p + σ2m.ρ[
          t].σ1p + σ2m.ρ[t].σ2p)) - 
  1/2 sh ((ρ[t].σ1m.σ1p + ρ[
         t].σ1m.σ2p + ρ[
         t].σ2m.σ1p + ρ[
         t].σ2m.σ2p) + \
(σ1m.σ1p.ρ[t] + σ1m.σ2p.ρ[
             t] + σ2m.σ1p.ρ[
             t] + σ2m.σ2p.ρ[t]) - 
         2 (σ1p.ρ[t].σ1m + σ1p.ρ[
              t].σ2m + σ2p.ρ[
              t].σ1m + σ2p.ρ[t].σ2m)) - 
      I Λ21 (σ1p.σ2m.ρ[t] - ρ[
           t].σ1p.σ2m) - 
      I Λ21 (σ2p.σ1m.ρ[t] - ρ[
           t].σ2p.σ1m) , ρ[0] == init}, 
   slnn, {t, 0, 1.6}, WorkingPrecision -> 44];
(*plot*)
Plot[{Evaluate[a[2, 2][t] /. solsq], Evaluate[a[4, 4][t] /. solsq], 
  Evaluate[a[1, 1][t] /. solsq], Evaluate[a[3, 3][t] /. solsq], 
  Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] + 
   Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq]}, {t,
   0, 1.1}, PlotRange -> {-2, 6}, 
 PlotLegends -> {"ud", "dd", "uu", "du", "sum"}]

Now I get the result:

which cannot be correct, because we can prove that:

1,"uu","ud","du","dd" must be between 0 and 1

2, "sum" should be a constant, it is the trace of matrix ρ, and it can be easily proved by taking the trace of both sides of our ODE. when we take the trace, we have

0.*10^-43 I a[1, 4][t] + 0.*10^-43 I a[4, 1][t]

If we increase the WorkingPrecision in the first line of the code, we can make it closer to 0. However, the plot still shows a very bad behavior no matter how I change the WorkingPrecision.

I try to pick the value of solution:

Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] +Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq] /. t -> 1
Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] +Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq] /. t -> 1.5
Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] +Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq] /. t -> 2

They give:

{1.00000000000000000000}
{1.000000000}
no significant digits available to display

That's why the solution seems so badly behaved, but I can't fix it. Increasing the working precision seems not to help. Can anyone help me solve this numerical issue? Thanks!