I have been trying to vectorize my code for the past few days to speed up root finding in a grid. Sorry in advance if my code is badly formatted. I have quite a bit of setup for the problem. It really starts here where I solve a large system of equations and find an interpolated function trcfnint:

```
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`

is only the vectorized version). This is very quick and easy, though `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 I do that, I only take the upper half bladder. The initial values for a are chosen because they are in a spectral gap on the a axis, and this saves a lot of time in the calculations. Generally I want the ability to quickly solve for many values of a. I investigated various questions about this vectorization: Findroot with a precompiled function with parameters, Jacobian for parallelized FindRoot with multiple variables, FindRoot with vector functions were the most useful. I can publish my successful vectorizations, but they suffer from being slower than the undisclosed implementation. If anyone has any advice, it would be greatly appreciated.