# How can I accelerate repeated interpolation of a very large dataset?

I am solving a 4D time-dependent PDE using a discretization method. My query is not about solving this; because I am doing it in the following code. The code is OK and it works fine. My main problem occurs when I want to use the solution of this PDE (its interpolation solution1[]) in a loop over the grid of points, i.e., the last piece of the written code. This part is quite slow.

ClearAll["Global*"];
t1 = AbsoluteTime[];
(***************PARAMETER SETTINGS***************)
TT = 5.; m1 = 10; m2 = 10; m3 = 10; m4 = 10;
size = m1*m2*m3*m4;
Print["The size of the system of ODEs is = ", size];
roRr = -0.4; roRz = -0.15; roRy = -0.15; rorz = 0.5; rory = 0.5; royz \
= 0.25;
sigmaR = 0.01; sigmar = 0.08; sigmaz = 0.1; sigmay = 0.4;
gammar = 0.08; gammaz = 0.08;
r = 0.02; θ = -210.; a = 0.08; b = 0.1; a1 = 0.08; b1 = 0.1; e \
= 1.15; κ = 0.0001;
(***************DOMAIN DISCRETIZATION***************)
Rmin = 0.; Rmax = 1.; rmin = 0.; rmax = 1.; ymin = -6.; ymax = 0.; \
zmin = 0.; zmax =(*14e*)4;
Print[{Rmin, Rmax}, "\[Cross]", {rmin, rmax}, "\[Cross]" {ymin, ymax},
"\[Cross]", {zmin, zmax}];

e1 = 0.45;
nx = xgrid1 = Range[Rmin, Rmax, (Rmax - Rmin)/(m1 - 1)];
ny = ygrid1 = Range[rmin, rmax, (rmax - rmin)/(m2 - 1)];
nz = zgrid1 = Range[ymin, ymax, (ymax - ymin)/(m3 - 1)];
nw = wgrid1 = Range[zmin, zmax, (zmax - zmin)/(m4 - 1)];
origrid = Flatten[Outer[List, nx, ny, nz, nw], 3];

(***************FILLING SOME MATRICES***************)
Idx = SparseArray[{{i_, i_} -> 1.}, {m1, m1}, 0]; Idy =
SparseArray[{{i_, i_} -> 1.}, {m2, m2}, 0];
Idz = SparseArray[{{i_, i_} -> 1.}, {m3, m3}, 0]; Idw =
SparseArray[{{i_, i_} -> 1.}, {m4, m4}, 0];
DR = KroneckerProduct[(SparseArray@DiagonalMatrix@nx), Idy, Idz, Idw];
Dr = KroneckerProduct[Idx, (SparseArray@DiagonalMatrix@ny), Idz, Idw];
Dy = KroneckerProduct[Idx, Idy, (SparseArray@DiagonalMatrix@nz), Idw];
Dz = KroneckerProduct[Idx, Idy, Idz, (SparseArray@DiagonalMatrix@nw)];
Id = KroneckerProduct[Idx, Idy, Idz, Idw];

Dy2 = KroneckerProduct[Idx, Idy, (SparseArray@DiagonalMatrix@Exp[nz]),
Idw];
(***************FOR FIRST SPATIAL VARIABLE***************)
hh = Differences[nx]; ww =
Table[hh[[i + 1]]/hh[[i]], {i, 1, Length[hh] - 1}];
{Min[hh], c = 6*Max[hh]}
A1d1 = Join[
Table[(ww [[
i - 1]] (-3 c^2 + hh[[i - 1]]^2 (-5 + 2 ww [[i - 1]])))/(
3 c^2 hh[[i - 1]] (1 + ww [[i - 1]])), {i, 2,
m1 - 1}], { -(1/hh[[Length[hh]]]) + hh[[Length[hh]]]/c^2}];
A1d2 = Join[{-(1/hh[[1]]) + hh[[1]]/c^2},
Table[-((2  hh[[i - 1]] (-1 + ww[[i - 1]]))/(3 c^2)) + (-1 +
ww[[i - 1]])/( hh[[i - 1]] ww[[i - 1]]), {i, 2, m1 - 1}], {1/
hh[[Length[hh]]]}];
A1d3 = Join[{1/hh[[1]]},
Table[(3/ww[[i - 1]] + (hh[[i - 1]]^2 (-2 + 5 ww[[i - 1]]))/c^2)/(
3 hh[[i - 1]] (1 + ww[[i - 1]])), {i, 2, m1 - 1}]];
(*2nd*)
A3d1 = Join[
Table[(3 hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) -
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (1 + (-3 + ww[[i - 1]]) ww[[i - 1]]))/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
hh[[i - 2]]/hh[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m1 - 1}], {0}];
A2d1 = Join[{(2 (3/hh[[1]]^2 + (5 + 2 (-2 + ww[[1]]) ww[[1]])/c^2))/(
3 (1 + ww[[1]]))}, Table[
((1 + hh[[i - 2]]/hh[[i - 1]]) (2 c^2 + 3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) - (2 c^2 +
3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
ww[[i - 1]]))
, {i, 3, m1 - 1}], {-4/c^2}];
A2d2 = Join[{-4/c^2, (
2 (-(3/hh[[1]]^2) + (-2 + ww[[1]] - 2 ww[[1]]^2)/c^2))/(
3 ww[[1]])}, Table[
1/(c^2 hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]]) ww[[
i - 1]]) (hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) - (1 + hh[[i - 2]]/
hh[[i - 1]]) (2 c^2 + hh[[i - 1]]^2 +
hh[[i - 1]]^2 (-1 + ww[[i - 1]]) ww[[i - 1]]))
, {i, 3, m1 - 1}], {2/c^2}];
A2d3 = Join[{2/c^2, (
6 c^2 + 2 hh[[1]]^2 (2 + ww[[1]] (-4 + 5 ww[[1]])))/(
3 c^2 hh[[1]]^2 ww[[1]] (1 + ww[[1]]))}, Table[
((1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (2 c^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) -
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
3 hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[i - 1]]^2 ww[[
i - 1]] (1 + ww[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m1 - 1}]];

dudx = Chop@
SparseArray[{Band[{1, 1}] -> A1d2, Band[{2, 1}] -> A1d1,
Band[{1, 2}] -> A1d3}, {m1, m1}];
d2udx2 = Chop@
SparseArray[{Band[{1, 1}] -> A2d2, Band[{2, 1}] -> A2d1,
Band[{1, 2}] -> A2d3, Band[{3, 1}] -> A3d1}, {m1, m1}];

(***************FOR SECOND SPATIAL VARIABLE***************)
hh = Differences[ny]; ww =
Table[hh[[i + 1]]/hh[[i]], {i, 1, Length[hh] - 1}];
{Min[hh], c = 6*Max[hh]}
A1d1 = Join[
Table[(ww [[
i - 1]] (-3 c^2 + hh[[i - 1]]^2 (-5 + 2 ww [[i - 1]])))/(
3 c^2 hh[[i - 1]] (1 + ww [[i - 1]])), {i, 2,
m2 - 1}], { -(1/hh[[Length[hh]]]) + hh[[Length[hh]]]/c^2}];
A1d2 = Join[{-(1/hh[[1]]) + hh[[1]]/c^2},
Table[-((2  hh[[i - 1]] (-1 + ww[[i - 1]]))/(3 c^2)) + (-1 +
ww[[i - 1]])/( hh[[i - 1]] ww[[i - 1]]), {i, 2, m2 - 1}], {1/
hh[[Length[hh]]]}];
A1d3 = Join[{1/hh[[1]]},
Table[(3/ww[[i - 1]] + (hh[[i - 1]]^2 (-2 + 5 ww[[i - 1]]))/c^2)/(
3 hh[[i - 1]] (1 + ww[[i - 1]])), {i, 2, m2 - 1}]];
(*2nd*)
A3d1 = Join[
Table[(3 hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) -
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (1 + (-3 + ww[[i - 1]]) ww[[i - 1]]))/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
hh[[i - 2]]/hh[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m2 - 1}], {0}];
A2d1 = Join[{(2 (3/hh[[1]]^2 + (5 + 2 (-2 + ww[[1]]) ww[[1]])/c^2))/(
3 (1 + ww[[1]]))}, Table[
((1 + hh[[i - 2]]/hh[[i - 1]]) (2 c^2 + 3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) - (2 c^2 +
3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
ww[[i - 1]]))
, {i, 3, m2 - 1}], {-4/c^2}];
A2d2 = Join[{-4/c^2, (
2 (-(3/hh[[1]]^2) + (-2 + ww[[1]] - 2 ww[[1]]^2)/c^2))/(
3 ww[[1]])}, Table[
1/(c^2 hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]]) ww[[
i - 1]]) (hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) - (1 + hh[[i - 2]]/
hh[[i - 1]]) (2 c^2 + hh[[i - 1]]^2 +
hh[[i - 1]]^2 (-1 + ww[[i - 1]]) ww[[i - 1]]))
, {i, 3, m2 - 1}], {2/c^2}];
A2d3 = Join[{2/c^2, (
6 c^2 + 2 hh[[1]]^2 (2 + ww[[1]] (-4 + 5 ww[[1]])))/(
3 c^2 hh[[1]]^2 ww[[1]] (1 + ww[[1]]))}, Table[
((1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (2 c^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) -
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
3 hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[i - 1]]^2 ww[[
i - 1]] (1 + ww[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m2 - 1}]];

dudy = Chop@
SparseArray[{Band[{1, 1}] -> A1d2, Band[{2, 1}] -> A1d1,
Band[{1, 2}] -> A1d3}, {m2, m2}];
d2udy2 = Chop@
SparseArray[{Band[{1, 1}] -> A2d2, Band[{2, 1}] -> A2d1,
Band[{1, 2}] -> A2d3, Band[{3, 1}] -> A3d1}, {m2, m2}];

(***************FOR THIRD SPATIAL VARIABLE***************)
hh = Differences[nz]; ww =
Table[hh[[i + 1]]/hh[[i]], {i, 1, Length[hh] - 1}];
{Min[hh], c = 6*Max[hh]}
A1d1 = Join[
Table[(ww [[
i - 1]] (-3 c^2 + hh[[i - 1]]^2 (-5 + 2 ww [[i - 1]])))/(
3 c^2 hh[[i - 1]] (1 + ww [[i - 1]])), {i, 2,
m3 - 1}], { -(1/hh[[Length[hh]]]) + hh[[Length[hh]]]/c^2}];
A1d2 = Join[{-(1/hh[[1]]) + hh[[1]]/c^2},
Table[-((2  hh[[i - 1]] (-1 + ww[[i - 1]]))/(3 c^2)) + (-1 +
ww[[i - 1]])/( hh[[i - 1]] ww[[i - 1]]), {i, 2, m3 - 1}], {1/
hh[[Length[hh]]]}];
A1d3 = Join[{1/hh[[1]]},
Table[(3/ww[[i - 1]] + (hh[[i - 1]]^2 (-2 + 5 ww[[i - 1]]))/c^2)/(
3 hh[[i - 1]] (1 + ww[[i - 1]])), {i, 2, m3 - 1}]];
(*2nd*)
A3d1 = Join[
Table[(3 hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) -
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (1 + (-3 + ww[[i - 1]]) ww[[i - 1]]))/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
hh[[i - 2]]/hh[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m3 - 1}], {0}];
A2d1 = Join[{(2 (3/hh[[1]]^2 + (5 + 2 (-2 + ww[[1]]) ww[[1]])/c^2))/(
3 (1 + ww[[1]]))}, Table[
((1 + hh[[i - 2]]/hh[[i - 1]]) (2 c^2 + 3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) - (2 c^2 +
3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
ww[[i - 1]]))
, {i, 3, m3 - 1}], {-4/c^2}];
A2d2 = Join[{-4/c^2, (
2 (-(3/hh[[1]]^2) + (-2 + ww[[1]] - 2 ww[[1]]^2)/c^2))/(
3 ww[[1]])}, Table[
1/(c^2 hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]]) ww[[
i - 1]]) (hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) - (1 + hh[[i - 2]]/
hh[[i - 1]]) (2 c^2 + hh[[i - 1]]^2 +
hh[[i - 1]]^2 (-1 + ww[[i - 1]]) ww[[i - 1]]))
, {i, 3, m3 - 1}], {2/c^2}];
A2d3 = Join[{2/c^2, (
6 c^2 + 2 hh[[1]]^2 (2 + ww[[1]] (-4 + 5 ww[[1]])))/(
3 c^2 hh[[1]]^2 ww[[1]] (1 + ww[[1]]))}, Table[
((1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (2 c^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) -
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
3 hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[i - 1]]^2 ww[[
i - 1]] (1 + ww[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m3 - 1}]];

dudz = Chop@
SparseArray[{Band[{1, 1}] -> A1d2, Band[{2, 1}] -> A1d1,
Band[{1, 2}] -> A1d3}, {m3, m3}];
d2udz2 = Chop@
SparseArray[{Band[{1, 1}] -> A2d2, Band[{2, 1}] -> A2d1,
Band[{1, 2}] -> A2d3, Band[{3, 1}] -> A3d1}, {m3, m3}];

(***************FOR FOURTH SPATIAL VARIABLE***************)
hh = Differences[nw]; ww =
Table[hh[[i + 1]]/hh[[i]], {i, 1, Length[hh] - 1}];
{Min[hh], c = 6*Max[hh]}
A1d1 = Join[
Table[(ww [[
i - 1]] (-3 c^2 + hh[[i - 1]]^2 (-5 + 2 ww [[i - 1]])))/(
3 c^2 hh[[i - 1]] (1 + ww [[i - 1]])), {i, 2,
m4 - 1}], { -(1/hh[[Length[hh]]]) + hh[[Length[hh]]]/c^2}];
A1d2 = Join[{-(1/hh[[1]]) + hh[[1]]/c^2},
Table[-((2  hh[[i - 1]] (-1 + ww[[i - 1]]))/(3 c^2)) + (-1 +
ww[[i - 1]])/( hh[[i - 1]] ww[[i - 1]]), {i, 2, m4 - 1}], {1/
hh[[Length[hh]]]}];
A1d3 = Join[{1/hh[[1]]},
Table[(3/ww[[i - 1]] + (hh[[i - 1]]^2 (-2 + 5 ww[[i - 1]]))/c^2)/(
3 hh[[i - 1]] (1 + ww[[i - 1]])), {i, 2, m4 - 1}]];
(*2nd*)
A3d1 = Join[
Table[(3 hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) -
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (1 + (-3 + ww[[i - 1]]) ww[[i - 1]]))/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
hh[[i - 2]]/hh[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m4 - 1}], {0}];
A2d1 = Join[{(2 (3/hh[[1]]^2 + (5 + 2 (-2 + ww[[1]]) ww[[1]])/c^2))/(
3 (1 + ww[[1]]))}, Table[
((1 + hh[[i - 2]]/hh[[i - 1]]) (2 c^2 + 3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) - (2 c^2 +
3 hh[[i - 1]]^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[
i - 1]]^2 (-1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (1 +
ww[[i - 1]]))
, {i, 3, m4 - 1}], {-4/c^2}];
A2d2 = Join[{-4/c^2, (
2 (-(3/hh[[1]]^2) + (-2 + ww[[1]] - 2 ww[[1]]^2)/c^2))/(
3 ww[[1]])}, Table[
1/(c^2 hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]]) ww[[
i - 1]]) (hh[[
i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])^2 (-1 +
ww[[i - 1]]) + (-1 + ww[[i - 1]]) (2 c^2 -
hh[[i - 1]]^2 ww[[i - 1]]) - (1 + hh[[i - 2]]/
hh[[i - 1]]) (2 c^2 + hh[[i - 1]]^2 +
hh[[i - 1]]^2 (-1 + ww[[i - 1]]) ww[[i - 1]]))
, {i, 3, m4 - 1}], {2/c^2}];
A2d3 = Join[{2/c^2, (
6 c^2 + 2 hh[[1]]^2 (2 + ww[[1]] (-4 + 5 ww[[1]])))/(
3 c^2 hh[[1]]^2 ww[[1]] (1 + ww[[1]]))}, Table[
((1 + (1 + hh[[i - 2]]/hh[[i - 1]])) (2 c^2 +
hh[[i - 1]]^2 (1 + hh[[i - 2]]/hh[[i - 1]])) -
hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/
hh[[i - 1]]) (3 + (1 + hh[[i - 2]]/hh[[i - 1]]))) ww[[
i - 1]] +
3 hh[[i - 1]]^2 (1 + (1 + hh[[i - 2]]/hh[[i - 1]])) ww[[
i - 1]]^2)/(c^2 hh[[i - 1]]^2 ww[[
i - 1]] (1 + ww[[i - 1]]) ((1 + hh[[i - 2]]/hh[[i - 1]]) +
ww[[i - 1]]))
, {i, 3, m4 - 1}]];

dudw = Chop@
SparseArray[{Band[{1, 1}] -> A1d2, Band[{2, 1}] -> A1d1,
Band[{1, 2}] -> A1d3}, {m4, m4}];
d2udw2 = Chop@
SparseArray[{Band[{1, 1}] -> A2d2, Band[{2, 1}] -> A2d1,
Band[{1, 2}] -> A2d3, Band[{3, 1}] -> A3d1}, {m4, m4}];

(****************BUILDING THE SYSTEM MATRIX********************)

B = SparseArray[
+(1/2 sigmaR^2 DR.(Id - DR)).KroneckerProduct[d2udx2, Idy, Idz, Idw]
+ (1/2 sigmar^2 Dr).KroneckerProduct[Idx, d2udy2, Idz, Idw]
+ (1/2 sigmay^2)*KroneckerProduct[Idx, Idy, d2udz2, Idw]
+ (1/2 sigmaz^2 Dz^2)*KroneckerProduct[Idx, Idy, Idz, d2udw2]
+ ((roRr*sigmaR*sigmar)*Sqrt[(DR.(Id - DR)).Dr]).KroneckerProduct[
dudx, dudy, Idz, Idw]
+ ((roRz*sigmaR*sigmaz)*(Dz.Sqrt[DR.(Id - DR)])).KroneckerProduct[
dudx, Idy, Idz, dudw]
+ ((rorz*sigmar*sigmaz)*(Dz.Sqrt[Dr])).KroneckerProduct[Idx, dudy,
Idz, dudw]
+ ((roRy*sigmaR*sigmay)*(Sqrt[DR.(Id - DR)])).KroneckerProduct[
dudx, Idy, dudz, Idw]
+ ((rory*sigmar*sigmay)*(Sqrt[Dr])).KroneckerProduct[Idx, dudy,
dudz, Idw]
+ ((royz*sigmay*sigmaz)*Dz).KroneckerProduct[Idx, Idy, dudz, dudw]
+ (a*(b*Id - DR)).KroneckerProduct[dudx, Idy, Idz, Idw]
+ (a1 (b1*Id - Dr)).KroneckerProduct[Idx, dudy, Idz, Idw]
+ ((r*Id - Dr).Dz).KroneckerProduct[Idx, Idy, Idz, dudw]
+ (κ (θ*Id - Dy)).KroneckerProduct[Idx, Idy, dudz,
Idw]
- (r*Id + Dy2)
- gammaz*(Dy2.(Dz.KroneckerProduct[Idx, Idy, Idz, dudw]))
];

payoff = Flatten@
Table[nw[[l]], {i, 1, m1}, {j, 1, m2}, {k, 1, m3}, {l, 1, m4}];

(*Imposing the boundaries*)
index1 = Flatten[
Table[{i, j, k, l}, {i, 1, m1}, {j, 1, m2}, {k, 1, m3}, {l, 1,
m4}], 3];

(*Boundary when R tends to Subscript[R, max]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m1}];
unk = D[Last[U12[t]], t];
unk2 = Drop[D[Take[U12[t], -4], t], -1];
Last[Simplify@
NDSolveFiniteDifferenceDerivative[2, nx, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == m1), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == m1 - 1), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete222 = DeleteCases[pos2, 0]; // AbsoluteTiming
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == m1 - 2), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == m1 - 3), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[3]] B[[delete222]] +
A121[[2]] B[[delete333]] +
A121[[1]] B[[delete444]]); // AbsoluteTiming

(*Boundary when R tends to Subscript[R, min]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m1}];
unk = D[First[U12[t]], t];
unk2 = Drop[D[Take[U12[t], 4], t], 1];
First[Simplify@
NDSolveFiniteDifferenceDerivative[2, nx, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == 1), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == 2), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete222 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == 3), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[1]] == 4), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[1]] B[[delete222]] +
A121[[2]] B[[delete333]] + A121[[3]] B[[delete444]]);

(*Boundary when r tends to Subscript[r, max]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m2}];
unk = D[Last[U12[t]], t];
unk2 = Drop[D[Take[U12[t], -4], t], -1];
Last[Simplify@
NDSolveFiniteDifferenceDerivative[2, ny, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == m2), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == m2 - 1), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete222 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == m2 - 2), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == m2 - 3), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[3]] B[[delete222]] +
A121[[2]] B[[delete333]] + A121[[1]] B[[delete444]]);

(*Boundary when r tends to Subscript[r, min]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m2}];
unk = D[First[U12[t]], t];
unk2 = Drop[D[Take[U12[t], 4], t], 1];
First[Simplify@
NDSolveFiniteDifferenceDerivative[2, ny, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == 1), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == 2), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete222 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == 3), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[2]] == 4), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[1]] B[[delete222]] +
A121[[2]] B[[delete333]] + A121[[3]] B[[delete444]]);

(*Boundary when y tends to Subscript[y, max]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m3}];
unk = D[Last[U12[t]], t];
unk2 = Drop[D[Take[U12[t], -4], t], -1];
Last[Simplify@
NDSolveFiniteDifferenceDerivative[2, nz, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A1210 = First[A111];
pos20 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == m3), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete1110 = DeleteCases[pos20, 0];
pos20 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == m3 - 1), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete2220 = DeleteCases[pos20, 0];
pos20 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == m3 - 2), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete3330 = DeleteCases[pos20, 0];
pos20 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == m3 - 3), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete4440 = DeleteCases[pos20, 0];
B[[delete1110]] = (A1210[[3]] B[[delete2220]] +
A1210[[2]] B[[delete3330]] + A1210[[1]] B[[delete4440]]);

(*Boundary when y tends to Subscript[y, min]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m3}];
unk = D[First[U12[t]], t];
unk2 = Drop[D[Take[U12[t], 4], t], 1];
First[Simplify@
NDSolveFiniteDifferenceDerivative[2, nz, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == 1), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == 2), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete222 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == 3), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[3]] == 4), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[1]] B[[delete222]] +
A121[[2]] B[[delete333]] + A121[[3]] B[[delete444]]);

(*Boundary when z tends to Subscript[z, max]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m4}];
unk = D[Last[U12[t]], t];
unk2 = Drop[D[Take[U12[t], -4], t], -1];
Last[Simplify@
NDSolveFiniteDifferenceDerivative[2, nw, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == m4), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == m4 - 1), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete222 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == m4 - 2), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == m4 - 3), var1[l] = l, var1[l] = 0], {l,
Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[3]] B[[delete222]] +
A121[[2]] B[[delete333]] + A121[[1]] B[[delete444]]);

(*Boundary when z tends to Subscript[z, min]*)
U12[t_] = Table[Subscript[u, i][t], {i, 1, m4}];
unk = D[First[U12[t]], t];
unk2 = Drop[D[Take[U12[t], 4], t], 1];
First[Simplify@
NDSolveFiniteDifferenceDerivative[2, nw, U12[t],
"DifferenceOrder" -> 2]];
% == 0;
FullSimplify@(unk /. Solve[D[%, t], unk]);
{b111, A111} = Normal@CoefficientArrays[%, unk2];
A121 = First[A111];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == 1), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete111 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == 2), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete222 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == 3), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete333 = DeleteCases[pos2, 0];
pos2 = Table[ind1 = index1[[l]];
If[(ind1[[4]] == 4), var1[l] = l, var1[l] = 0], {l, Length@index1}];
delete444 = DeleteCases[pos2, 0];
B[[delete111]] = (A121[[1]] B[[delete222]] +
A121[[2]] B[[delete333]] + A121[[3]] B[[delete444]]);
8

B = SparseArray[B];

k1 = .05;
Monitor[
lines = NDSolve[
{D[v[t], t] == B.v[t], v[0] == initc[[All, 2]]}, v[t], {t, 0, TT},
Method -> {"FixedStep", "StepSize" -> k1,
Method -> {"ExplicitRungeKutta",
"DifferenceOrder" -> 4, "StiffnessTest" -> False}},
PrecisionGoal -> 5, AccuracyGoal -> 5,
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]; // AbsoluteTiming

s = v[t] /. lines[[1]];
solution1[time_] := s /. t -> time
sol2 = solution1[TT];

ClearAll[Idx, Idy, Idz, Idw, dudx, dudy, dudz, d2udx2, d2udy2, d2udz2,
hh, ww, delete, delete1, delete11];
set1 = Chop@
Flatten[Table[{Flatten@{origrid[[i]], sol2[[i]]}}, {i, 1,
Length[payoff]}], 1];
gg = Interpolation@set1;

NN = 10;
m = 120;
tInterval = Range[0, TT, TT/m];
h = (TT/m)/NN;

leng = Length[payoff];
tab1 = Table[i, {i, size}];

Aa = Parallelize@Table[
setsol1 = Table[
g1 = With[{current = solution1[input]},
Interpolation@Threshold@Flatten[

Table[{Flatten@{origrid[[i]], (current)[[i]]}}, {i,
leng}], 1]
];
g1[e1, 0.03, -4.089, e]
,
{input, Table[(tInterval[[j]] + k*h), {k, NN}]}
];
h*Sum[setsol1[[k]], {k, NN}]
, {j, m}]; // AbsoluteTiming

t2 = AbsoluteTime[] - t1;
Print["The whole computational time is = ", t2];


I know this code is large but my problem is only about its last part (when $$m=120$$, i.e., for calculating the list for Aa).

How can we accelerate with Interpolation[] inside a loop? In other words, how can we make fast repeated interpolations on very large 4D grids?

I have now used the parallel command ParallelTable[] which requires 14.5 seconds on my system. Noting that in my actual tests, I wish to increase the number of discretization points of the PDE which lead to higher CPU time. Therefore, I will be thankful if someone could give me some hint to improve the performance. Can we "Compile[]" a list of vectors? Can we do some part in the GPU using Mathematica CUDA applicability?

• What do you calculate with this code, what should be the result? Oct 5, 2019 at 19:16
• After solving the PDE, I must find the values of the solution at a special location of the domain (e1, 0.03, -4.089, e) for different times. This must be done $m$ times which $m$ comes from a financial application. Here, it shows bi-monthly payment of a quanto CDS par spread. So, the final result of the Aa is a list of numbers/values/par-spreads, which I finally do some operations (such as addition) on this list.
– Faz
Oct 5, 2019 at 19:22
• I'm voting to close this question as off-topic because OP does not care at all to produce a minimal example. Oct 5, 2019 at 23:27
• @HenrikSchumacher I think you can optimize this problem. Oct 6, 2019 at 11:25
• @Fazlollah Can you explain the origin of this problem? Are there any published articles on this topic? Oct 6, 2019 at 11:29

This is the actual problem, condensed in only a couple of lines:

m1 = 10; m2 = 10; m3 = 10; m4 = 10;
dims = {10, 10, 10, 10};
grids = {
Subdivide[0., 1., m1 - 1],
Subdivide[0., 1., m2 - 1],
Subdivide[-6., 0., m3 - 1],
Subdivide[0., 4., m4 - 1]
};
origrid = Tuples[grids];
leng = Length[origrid];
pt = DeveloperToPackedArray[{0.45, 0.03, -4.089, 1.15}];

vals = RandomReal[{-1, 1}, leng];
g = Interpolation@ Flatten[Table[{Flatten@{origrid[[i]], vals[[i]]}}, {i, leng}], 1]; // AbsoluteTiming // First
result = g @@ pt; // AbsoluteTiming // First


0.04719

0.000149

That is: We need to evaluate a certain InterpolatingFunction g computed from values vals on a regular, rectangular grid origrid only at a single fixed point pt. The values of vals may change over time, but the point pt is fixed. As one can see, creating the interpolating function takes about 300 times as long as evaluating it.

A first step is to replace the suspiciously complicated way of building the array that is fed to Interpolation:

g = Interpolation@Join[origrid, Partition[vals, 1], 2]; // AbsoluteTiming // First


0.011582

This speeds things up already by a factor of 4.

# Avoiding Interpolation

Denote the the grid points by $$P_{i_1,i_2,i_3,i_4}$$ with $$1 \leq i_k \leq n_k$$ and the functions values that we want to interpolate by $$v_{i_1,i_2,i_3,i_4}$$ with $$1 \leq i_k \leq n_k$$. Denote by $$x$$ a point in $$\mathbb{R}^4$$ on which we want to evaluate the function $$g$$, where the latter satisfies $$g(P_{i_1,i_2,i_3,i_4}) = v_{i_1,i_2,i_3,i_4}$$.

It depends somehow on the interpolation scheme that one wants to use. IIRC, Interpolation uses tensor product splines of order $$3$$ by default. That means, the stencil of a point evalutation has size $$4 = 3 + 1$$ per dimension. So what one should do:

The point $$x$$ lies in a certain cell of the rectangular grid. Denote the grid point to the "lower left" by $$P_{j_1,j_2,j_3,j_4}$$.

Then $$g(x)$$ depends only on the values $$v_{i_1,i_2,i_3,i_4}\quad \text{with} \quad j_k -1 \leq i_k \leq j_k+2, \quad k \in \{1,\dotsc,4\}.$$

More precisely, one has $$g(x) = \sum_{j_k -1 \leq i_k \leq j_k+2} \omega_{i_1,i_2,i_3,i_4} \, v_{i_1,i_2,i_3,i_4}, \quad \text{where} \quad \omega_{i_1,i_2,i_3,i_4} := \varphi_{i_1,i_2,i_3,i_4}(x)$$ and where $$\varphi_{i_1,i_2,i_3,i_4}$$ is the interpolating function that satisfies $$\varphi_{i_1,i_2,i_3,i_4}(P_{i_1,i_2,i_3,i_4}) = 1$$ and that vanishes on all other grid points. The weights $$\omega_{i_1,i_2,i_3,i_4}$$ have to be calculcated once and only for the grid points in the stencil, hence on only $$256= (3+1)^4$$ points.

Since tensor product splines are employed, one has $$\omega_{i_1,i_2,i_3,i_4} = \varphi_{i_1,i_2,i_3,i_4}(x) = \varphi_{1,i_1}(x_1) \, \varphi_{2,i_2}(x_2) \, \varphi_{3,i_3}(x_3) \, \varphi_{4,i_4}(x_4)$$ where $$\varphi_{k,i}$$ is the $$i$$-th interpolating in the $$k$$-th direction.

This can be implemented as follows:

w = Flatten[
SparseArray@Interpolation[Transpose[{#1, IdentityMatrix[#2, WorkingPrecision -> MachinePrecision]}]]@#3 &,
{grids, dims, pt}
]
];
idx = Flatten[w["NonzeroPositions"]];
w = w["NonzeroValues"];


Now, $$g(x)$$ can be evaluated as follows:

result2 = vals[[idx]].w; // AbsoluteTiming // First
result2 - result


0.000015

0.

Thus, the evaluatation is about 3000 times faster than generating the InterpolatingFunction g even 10 times faster than just evaluating it.

It becomes even more impressive if one performs evaluations in bulk, so that machine acceleration of matrix-vector products can be exploited (and Mathematica's interpretation overhead can be avoided):

vals = RandomReal[{-1, 1}, {120, leng}];
result = Table[
With[{current = vals[[j]]},
Interpolation[Flatten[Table[{Flatten@{origrid[[i]], current[[i]]}}, {i, leng}], 1]] @@ pt
],
{j, 1, 120}
]; // AbsoluteTiming // First
result2 = vals[[All, idx]].w; // AbsoluteTiming // First
Max[Abs[result - result2]]


5.25747

0.000176

4.996*10^-16

This is almost a 30000-fold speed-up.

After redefining solution1 appropriately, list Aa can then be obtained this way:

s = NDSolveValue[
{
D[v[t], t] == B.v[t],
v[0] == initc[[All, 2]]
},
v,
{t, 0, TT},
Method -> {
"FixedStep", "StepSize" -> k1,
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"StiffnessTest" -> False}
},
PrecisionGoal -> 5,
AccuracyGoal -> 5
];
(* Creating an interpolation function that returns only the required entries of s *)
solution1 = ReplacePart[s, 4 -> s[[4, All, All, idx]]];
Aafast = Dot[
solution1[Outer[Plus, tInterval[[1 ;; m]], Range[h, NN h, h]]],
w,
ConstantArray[h, NN]
]; // AbsoluteTiming // First
Max[Abs[Aafast - Aa]]


0.009964

1.38778*10^-17

On my machine, this runs about 2200 times faster than the original code for Aa. The new definition of solution1 enforces that only the relevant entries of the large vector a computed. Anyways, solution1 is still the bottleneck here.

• Fantastic 5436x acceleration on my laptop. Henrik, you are a genius! Oct 6, 2019 at 21:15
• A step by step useful optimization of the code. This is indeed great Henrik. Thanks a lot. I got a query. Why the length of each output/vector coming from your solution1[] is 256 (even when we change $m_1,m_2,m_3,m_4$)! Consider that I wish to plot the solution at the final time solution1[TT], so in my old and slow implementation, the length of the code was as same as the length of the origrid (for any $m_i$) and I was able to plot the solution while now I cannot! Can you please write some hints about this as well.
– Faz
Oct 7, 2019 at 6:02
• +1. So, you voted to close and then answered anyways ;) That's what I call an addict... perhaps we need a Alco-SE-lics Anonymous ;-) Well done. Oct 7, 2019 at 10:00
• @user21 It was a long, boring weekend you know... ;) Oct 7, 2019 at 12:09
• @Fazlollah The 256 comes from the fact that the stencil of the interpolation has size $4 \times 4 \times 4 \times 4$. That depends only on the order of the iterpolation scheme, not on the size of the grid. But I changed the code so that the result of NDSolveValue is stored in s (this is now what the old solution1 was). So if you want to plot the solution, just use s instead. Oct 7, 2019 at 12:12