I want to evaluate integrals of the kind
$$ \int_0^1 x \tanh{x}\, \mathcal{J}_{l + k}(\alpha_{n_1,l+k}\,x)\mathcal{J}_{l}(\alpha_{n_2,l}\,x) \, dx $$ where $\alpha_{n,l}$ is the $n^{th}$ zero of the $l^{th}$ Bessel function, $\mathcal{J}_l(x)$.
Currently, I'm using
Integ[l_, k_, n1_, n2_] := NIntegrate[x Tanh[x] BesselJ[l + k, x bjz[l + k, n1]] BesselJ[l, x bjz[l, n2]], {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}, WorkingPrecision -> 6, PrecisionGoal -> 5, AccuracyGoal -> 10];
where I define
bjz[n_?NumericQ, k_?NumericQ] := bjz[n, k] = N[BesselJZero[n, k]]
Let's say $k = 1$ here. For small n's and $l$, this method works fast enough for my purposes but say I take $l = 100, n_1 = 100, n_2 = 100$. Then it takes around 1 sec.
I tried using "LevinRule" and while that does speed up the integrals at large $l$ and $n$, it's slower than the other method at smaller values.
What would be the fastest way of evaluating such integrals such that it works for $0\leq l \leq 100$ and $0\leq n \leq 100$? I'm primarily interested in $k = 1,2$, although that shouldn't change much. There must be a general method for evaluating integrals of Bessel functions over a finite range, since these are pretty ubiquitous.