How to speed up (or parallelize) NDSolve?

Question

Is there a way to speed up the following code:

m = 1;
w = (a^2 + (m/p)^2)^(1/2); 
z = (b^2 - a^2)/w^2; 
y[z_] = ArcTan[z^(1/2)]/z^(1/2);
R[2, 4, 0, b_?NumericQ, a_?NumericQ] = (3 + 2 z - 3 (1 + z) y[z])/(z^2 w^3);
R[2, 2, 1, b_?NumericQ, a_?NumericQ] = (-3 + (3 + z) y[z])/(z^2 w^3);

II[n_, q_, r_, b_?NumericQ, a_?NumericQ] := ((b^(2 q + 2)) (a^(r + 1)) )/(4 \[Pi]^2 (2 q)!!) NIntegrate[p^(n + 1) R[n, q, r, b, a] E^-(p^2 + m^2)^(1/2), {p, 0, \[Infinity]}];

t0 = 0.1; tend = 1;
AbsoluteTiming[
 s1 = NDSolve[{D[II[2, 4, 0, b[t], a[t]], t] == a[t], D[II[2, 2, 1, b[t], a[t]], t] == b[t], b[t0] == 1, a[t0] == .1}, 
               {b[t], a[t]}, {t, t0, tend}]
]

Output:

Is it possible to parallelize NDSolve since that's where most of the time is spent?

This is a demo code. The actual code takes much longer (~30 mins).

Answer 1

We can reduce computational time by 25 times with extended expressions for coefficients of equations and by reducing integration time with option PrecisionGoal -> 4 as follows

Clear["Global`*"]

m = 1;
w = (a^2 + (m/p)^2)^(1/2);
z = (b^2 - a^2)/w^2;
y1 = ArcTan[z^(1/2)]/z^(1/2);
(*R[2,4,0,b_,a_]:=*)R1 = (3 + 2 z - 3 (1 + z) y1)/(z^2 w^3);
(*R[2,2,1,b_,a_]:=*)R2 = (-3 + (3 + z) y1)/(z^2 w^3);
(*II[n_,q_,r_,b_?NumericQ,a_?NumericQ]*)

A1 = (((b^(2 q + 2)) (a^(r + 1)))/(4 \[Pi]^2 (2 q)!!) p^(n + 
       1) R1 E^-(p^2 + m^2)^(1/2)) /. {n -> 2, q -> 4, r -> 0}; A11 = 
 D[A1, a]; A12 = D[A1, b];

a11[x_?NumericQ, y_?NumericQ] := 
  NIntegrate[
   If[x != y, A11 /. {a -> x, b -> y}, (
    E^-Sqrt[1 + p^2] p^5 x^10 (7 - 8 p^2 x^2))/(
    26880 Sqrt[1/p^2 + x^2] (\[Pi] + p^2 \[Pi] x^2)^2)], {p, 10^-8, 
    40}, PrecisionGoal -> 4];

a12[x_?NumericQ, y_?NumericQ] := 
  NIntegrate[
   If[x != y, A12 /. {a -> x, b -> y}, (
    E^-Sqrt[1 + p^2] p^7 x^10 Sqrt[1/p^2 + x^2] (35 + 32 p^2 x^2))/(
    13440 \[Pi]^2 (1 + p^2 x^2)^3)], {p, 10^-8, 40}, 
   PrecisionGoal -> 4];

A2 = (((b^(2 q + 2)) (a^(r + 1)))/(4 \[Pi]^2 (2 q)!!) p^(n + 
       1) R2 E^-(p^2 + m^2)^(1/2)) /. {n -> 2, q -> 2, r -> 1}; A21 = 
 D[A2, a]; A22 = D[A2, b];

a21[x_?NumericQ, y_?NumericQ] := 
 NIntegrate[
  If[x != y, A21 /. {a -> x, b -> y}, (
   E^-Sqrt[1 + p^2] p^5 x^7 (14 + 5 p^2 x^2))/(
   840 Sqrt[1/p^2 + x^2] (\[Pi] + p^2 \[Pi] x^2)^2)], {p, 10^-8, 40}, 
  PrecisionGoal -> 4]; 
a22[x_?NumericQ, y_?NumericQ] := 
 NIntegrate[
  If[x != y, A22 /. {a -> x, b -> y}, (
   E^-Sqrt[1 + p^2] p^7 x^7 Sqrt[1/p^2 + x^2] (7 + 5 p^2 x^2))/(
   140 \[Pi]^2 (1 + p^2 x^2)^3)], {p, 10^-8, 40}, PrecisionGoal -> 4];

eqs = {a11[x[t], y[t]] x'[t] + a12[x[t], y[t]] y'[t] - x[t] == 0, 
   a21[x[t], y[t]] x'[t] + a22[x[t], y[t]] y'[t] - y[t] == 0};
t0 = 0.1; tend = 1;
AbsoluteTiming[
 sol = NDSolve[{eqs, y[t0] == 1, x[t0] == .1}, {x[t], y[t]}, {t, t0, 
    tend}]]

Note, that we need exact expressions for a11,a12,a21,a22 as a limit at a->b. Computation time is

(*{14.3547, {{x[t] ->...}}}*)

Visualization

Plot[Evaluate[{x[t], y[t]} /. sol[[1]]], {t, t0, tend}, 
 AxesLabel -> Automatic, PlotLegends -> {"a", "b"}]