Repeated calling of a ParametricNDSolve result in a faster way?

Question

I recently learned of ParametricNDSolve and found that I need to compute it once and then call its result sol at different parameters x,y from rmat many times in the following example code. But always only evaluated at the same grid ttlist of the ODE time variable t, if this is relevant at all.

It turns out to be the bottleneck of my task. In reality I have much larger Ntt and Nr, so I need to speed it up as much as possible. (The form of eqs doesn't matter here.) Therefore I was wondering if any speed optimization possible? Either using Compile or not is fine.

Not sure if this is a stupid question as the result of ParametricNDSolve is probably not compilable. Just out of hope and curiosity...

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;
ttmax = tmax; ttmesh = ttmax/4; ttlist = Table[tt, {tt, -ttmax, ttmax, ttmesh}];
Ntt = Length@ttlist; Nr = 4;
rmat = RandomReal[1, {Nr, Nr, Ntt, Ntt, 2}];
sol = ParametricNDSolveValue[{eqs, ic}, Outer[#2[#] &, ttlist, var], {t, -tmax, tmax}, Element[#, Reals] & /@ {x, y}];
data = ParallelTable[
   With[{rxy = rmat[[ix, iy]]}, 
    Module[{Amat = 
       Table[Apply[sol, rxy[[it1, it2]]][[it1]], {it1, Ntt}, {it2, Ntt}]}, 
Total[Amat, 2](example of much faster computation that follows)]], {ix, Nr}, {iy, Nr}];
data // Abs // MinMax