1function z = wv_bjzeros(m, k)
2% WV_BJZEROS The first k positive zeros of the Bessel function J_m.
3%
4% Bracketed by a scan for sign changes and polished by Newton's method with
5% J_m' = (J_{m-1} - J_{m+1})/2, so each is good to a unit or two in the last
6% place. These are the exact frequencies of a circular drum, up to a factor
7% c/(2 pi R), and nothing about them is approximate.
8z = zeros(1, k);
9x = max(m, 0.5) : 0.05 : (m + 4 * k + 10);
10s = besselj(m, x);
11found = 0;
12for i = 1:numel(x) - 1
13 if s(i) * s(i+1) < 0 && found < k
14 r = x(i) - s(i) * (x(i+1) - x(i)) / (s(i+1) - s(i));
15 for it = 1:20
16 d = 0.5 * (besselj(m - 1, r) - besselj(m + 1, r));
17 st = besselj(m, r) / d;
18 r = r - st;
19 if abs(st) < 1e-15 * r
20 break;
21 end
22 end
23 found = found + 1;
24 z(found) = r;
25 end
26end
27end