During the process of evaluating the PSD from Sec. V of Liu and Thorn, I chanced to write a simpler root finding algorithm applying bisection to find the roots of J1(x).

The difference between this and the algorithm by

Greg von Winckel goes as

The difference is of the order 10^-7 and can be reduced by reducing the tolerance. However, It should however be noted that bisection is a crude algorithm for rough usage and differences become pronounced for larger n.

function X = bessel_zeros(k, n)
%BESSEL_ZEROS : calculates the first k zeros of the bessel function Jn
%   k : No. of ROOTS to be evaluated
%   n : Order of the Bessel function Jn
%   X : stores the roots in serial order i.e. X(1) gives the first root,
%   X(2) gives the second and so on
X = zeros(1,k);% empty array of lenght k
count = 1;%this acts as a counter to k
dx = 0.1;%step size within which, by assumption, no roots exist
x = 0;