How to speed up by Compile?

I need to NDSolve a system many times by scanning some of its parameters and do some matrix calculation with the (discretized) solutions. The example is the following code. The system eqs looks formidable, but the concrete form doesn’t matter.

Because the computation is quite heavy (I actually need much denser meshes for scanning the parameters, I mean the lists I defined in the middle), I was wondering if any further speed improvement possible. Probably by properly Compile the program? I’ve heard of great gain by compile to C or machine code. But I somehow got lost in the related documentation. Any suggestion would be appreciated.

Note that the NDSolve part sometimes takes less or comparable time than the matrix calculations followed, depending on how dense we scan particular parameters.

eqs = {I Derivative[1][B1][t] == 
    Abs[B1[t]]^2 (m B1[t] + B3[t] (x - I y + Cos[t] - I Sin[t])) + 
     Abs[B2[t]]^2 (m B1[t] + 
        B3[t] (x - I y + Cos[t] - I Sin[t])) + (B1[t] Conjugate[
          B3[t]] + B2[t] Conjugate[B4[t]]) (-m B3[t] + 
        B1[t] (x + I y + Cos[t] + I Sin[t])), 
   I Derivative[1][B2][t] == 
    Abs[B1[t]]^2 (m B2[t] + B4[t] (x - I y + Cos[t] - I Sin[t])) + 
     Abs[B2[t]]^2 (m B2[t] + 
        B4[t] (x - I y + Cos[t] - I Sin[t])) + (B1[t] Conjugate[
          B3[t]] + B2[t] Conjugate[B4[t]]) (-m B4[t] + 
        B2[t] (x + I y + Cos[t] + I Sin[t])), 
   I Derivative[1][B3][t] == 
    B1[t] (m (B3[t] Conjugate[B1[t]] + 
           B4[t] Conjugate[B2[t]]) + (Abs[B3[t]]^2 + 
           Abs[B4[t]]^2) (x + I y + Cos[t] + I Sin[t])) + 
     B3[t] (-m (Abs[B3[t]]^2 + 
           Abs[B4[t]]^2) + (B3[t] Conjugate[B1[t]] + 
           B4[t] Conjugate[B2[t]]) (x + Cos[t] - I (y + Sin[t]))), 
   I Derivative[1][B4][t] == 
    B2[t] (m (B3[t] Conjugate[B1[t]] + 
           B4[t] Conjugate[B2[t]]) + (Abs[B3[t]]^2 + 
           Abs[B4[t]]^2) (x + I y + Cos[t] + I Sin[t])) + 
     B4[t] (-m (Abs[B3[t]]^2 + 
           Abs[B4[t]]^2) + (B3[t] Conjugate[B1[t]] + 
           B4[t] Conjugate[B2[t]]) (x + Cos[t] - I (y + Sin[t])))};
ic = {B1[-tmax] == 1, B2[-tmax] == 0, B3[-tmax] == 0, B4[-tmax] == 1};
var = {B1, B2, B3, B4};

tmax = 25; m = 0.4;
tlist = Table[i, {i, 10, 30, 5}];
ttmax = tmax; ttmesh = ttmax/10; ttlist = Table[tt, {tt, -ttmax, ttmax, ttmesh}];
rmax = 2.0; rmesh = 0.1; rlist = Table[r, {r, -rmax, rmax, rmesh}];
omax = 1.0; omesh = omax/60; olist = Table[o, {o, -omax, omax, omesh}];
smat = Table[Exp[-(t1 - myt)^2/5 + I o t1], {o, olist}, {myt, tlist}, {t1, ttlist}];
G0 = {{x^2, y}, {-y, x y}};

data = ParallelTable[
   Module[{sol = Last[#] & /@ (First@NDSolve[{eqs, ic}, var, {t, -tmax, tmax}]), Bsol},
    Bsol = Table[Map[#@t1 &, sol, {1}], {t1, ttlist}]; 
    Map[#.G0.ConjugateTranspose[#] &@(ttmesh Partition[#1.Bsol, 2]) &, smat, {2}]], {x, rlist}, {y, rlist}];

tst = 
   Block[{x = rlist[[1]], y = rlist[[1]]}, 
    Module[{sol = (Last[#1] &) /@ First[NDSolve[{eqs, ic}, var, {t, -tmax, tmax}]], 
      Bsol}, Bsol = Table[Map[#1[t1] &, sol, {1}], {t1, ttlist}]; 
     Map[(#1 . G0 . ConjugateTranspose[#1] &)[ttmesh Partition[#1 . Bsol, 2]] &, 
      smat, {2}]]]; // AbsoluteTiming
(* {0.0249631, Null} *)

pnd = ParametricNDSolveValue[{eqs, ic}, 
   ArrayReshape[Outer[#2[#] &, ttlist, var], {Length@ttlist, 2, 2}], {t, -tmax, 
    tmax}, {x, y}];

cmap = With[{ttmesh = ttmesh, G0 = G0, smat = smat}, 
   Compile[{{x, _Real}, {y, _Real}, {Bsol, _Complex, 3}(*,{smat,_Complex,3}*)}, 
    With[{mat = ttmesh smat . Bsol}, 
     Module[{m, n}, {m, n} = Dimensions[mat][[;; 2]]; 
      Table[(mat[[i, j]] . G0 . ConjugateTranspose[mat[[i, j]]]), {i, m}, {j, n}]]], 
    CompilationTarget -> C]];

tst3 = Block[{x = rlist[[1]], y = rlist[[1]]}, 
    With[{Bsol = pnd[x, y]}, cmap[x, y, Bsol(*,smat*)]]]; // AbsoluteTiming
(* {0.0196469, Null} *)

tst3 - tst // Abs // Max
(* 2.8917*10^-14 *)