This paper describes an algorithm to calculate a large number of roots of the cross-product of Bessel functions and of their first derivatives. The algorithm initially finds the roots of the zeroth order using an auxiliary function that exhibits the same roots as the original cross-products but with better behavior for numerical root search with the Newton-Raphson algorithm. In order to find the roots for higher orders, the algorithm follows a pyramidal scheme using the interlacing property of the cross-product of Bessel functions. The algorithm shows globally convergent behavior for a large range of values of the argument and of the order of the Bessel functions. The roots can be computed to any precision, limited only by the computer implementation, and the convergence is attained in six iterations per root in average, showing a much better performance than previous works for the calculation of these roots.