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