# Repeated calling of a ParametricNDSolve result in a faster way?

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

• As I previously advised, if you're interested in compiling stuff, you should familiarize yourself with this. Currently, ParametricFunction[] (which is the head of what ParametricNDSolveValue[] returns) isn't in there. Feb 18 at 13:08
• @J.M. Yes, I note that in my question. But still hoping for any possible improvement. Feb 18 at 13:10
• @xiaohuamao It is not slow code, it takes seconds on my laptop. How you supposed improve it, to milliseconds? Feb 19 at 16:03
• @xiaohuamao Can you tell what actually do you try to compute? Your data looks like {{1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}}. How we can compare faster code with this one? Feb 20 at 23:35
• @AlexTrounev Sorry, the bottleneck part is the calculation of Amat inside Module inside data. I need this and then some other computation that is much faster and hence I just put 1 as a placeholder. Feb 21 at 2:30