I spent the last few days trying to vectorize my code to speed up root search on a grid. Sorry in advance if my code is badly formatted. I have a good amount of configuration for the problem. It really starts from here, where I solve a great system of equations and find an interpolated trcfnint function:

```
Clear("`*");
E2H(E_, c_) := 1/(c^2 - 1) (E)
E2H1 = Compile({{E, _Real}, {c, _Real}}, 1/(c^2 - 1) (E))
c = 1.1
E1 = 0.5
H = E2H1(E1, c)
If(c < 1 ,
period = Abs(
2 NIntegrate(
1/Sqrt(2 H - 1/(c^2 - 1) (u^4/2 - u^2)), {u, -Sqrt(
1 - Sqrt(1 + 4 (c^2 - 1) H)),
Sqrt(1 - Sqrt(1 + 4 (c^2 - 1) H))})),
period = Abs(
2 NIntegrate(
1/Sqrt(2 H - 1/(c^2 - 1) (u^4/2 - u^2)), {u, -Sqrt(
1 + Sqrt(1 + 4 (c^2 - 1) H)),
Sqrt(1 + Sqrt(1 + 4 (c^2 - 1) H))})));
sgtwv = NDSolveValue({(c^2 - 1) u2''(t) + u2(t)^3 - u2(t) == 0,
u2(0) == 0, u2'(0) == Sqrt(2 H) },
u2, {t, 0, 50});
grid = Table(a + I b, {a, -1, 1, 0.02}, {b, -1, 1, 0.02});
gvals = Riffle(Flatten(grid), Flatten(grid));
A(t_) := SparseArray(
Flatten(Table({{2 j - 1, 2 j} ->
1, {2 j,
2 j - 1} -> -1/(c^2 - 1) (gvals((j))^2/(c^2 - 1) +
3 sgtwv(t)^2 - 1)}, {j, Length(gvals)}), 1))
initialc =
Table(If(Mod(i, 4) == 1 || Mod(i, 4) == 0, 1, 0), {i,
2 Length(gvals)});
ff = NDSolveValue({q'(t) == A(t).q(t), q(0) == initialc}, {q}, {t, 0,
3 period}); // AbsoluteTiming
ggg = Flatten(Through(ff(period)));
trc = ggg((1 ;; -1 ;; 4)) + ggg((4 ;; -1 ;; 4));
listint = ArrayReshape(trc, Dimensions(grid));
trcfnint = ListInterpolation(listint, {{-1, 1}, {-1, 1}})
evs1(a_?VectorQ, b_?VectorQ, d_?VectorQ) :=
trcfnint(a, b) - 2 Cos(c period (a + I b)/(c^2 - 1) + d);
evs2(a_, b_, d_) :=
trcfnint(a, b) - 2 Cos(c period (a + I b)/(c^2 - 1) + d);
```

My goal is to find (a, b, d) in the interpolated range so that `evs2`

is 0 (`evs1`

it's just the vectorized version). This is very fast and easy when `b = 0`

:

```
imAxis1 = {}; Do(
plot1 = Plot(
trcfnint(a, 0) - 2 Cos(c period a/(c^2 - 1) + (Theta)), {a, 0, 1},
Mesh -> {{0}}, MeshFunctions -> {#2 &},
MeshStyle -> PointSize(Medium), PlotRange -> {{0, 1}, {0, 0}},
PlotPoints -> 100);
AppendTo(imAxis1,
Sort@Cases(Normal@plot1, Point({x_, y_}) -> x,
Infinity)), {(Theta), 0, 2 Pi + 0.1, 0.01});
imAxis1 = Flatten(imAxis1); imAxis1 =
Transpose({ConstantArray(0, Length(imAxis1)), imAxis1});
ListPlot(imAxis1, PlotStyle -> Directive(Red, PointSize(0.005)),
PlotRange -> {{-1, 1}, {-1, 1}},
FrameLabel -> {Style("Re((Lambda))", FontSize -> 20),
Style("Im((Lambda))", FontSize -> 20)},
LabelStyle -> Directive(Black), ImageSize -> Large, AspectRatio -> 1,
Frame -> True, Axes -> False)
```

When `b != 0`

, the most reliable results I have are using the non-vectorized function `evs2`

in the following code:

```
total = {};
Table(Do(
root1 =
FindRoot({Re(evs2(aa, b, dd)) == 0,
Im(evs2(aa, b, dd)) == 0}, {{aa, a}, {dd, theta}},
MaxIterations -> 20, AccuracyGoal -> 6);
test = Abs(evs2(aa, b, dd)) /. root1;
If(Abs(test) < 10^(-5), AppendTo(total, {b, root1((1))((2))}));
Break, {theta, 0, 2 Pi, 2}), {a, {0.5, 0.6}}, {b, 0, 0.18,
0.001}); // AbsoluteTiming
ListPlot(total, PlotStyle -> Directive(Red, PointSize(0.005)),
PlotRange -> {{-1, 1}, {-1, 1}},
FrameLabel -> {Style("Re((Lambda))", FontSize -> 20),
Style("Im((Lambda))", FontSize -> 20)},
LabelStyle -> Directive(Black), ImageSize -> Large, AspectRatio -> 1,
Frame -> True, Axes -> False)
```

When executing this, I only take the top half bubble. The initial values for a are chosen because they are in a spectral space on the axis a and this saves a lot of time in the calculations. In general, I would like to have the ability to quickly solve many values of a. I've examined several different questions about this vectorization: Findroot with a precompiled function with parameters, Jacobian for FindRoot in parallel with multiple variables, FindRoot with vector functions were the most useful. I can publish my successful vectorizations, however they suffer being slower than the non-vectorized implementation. If anyone has any advice I would greatly appreciate it.