function rj0 = jn_zeros ( n, nt )

%*****************************************************************************80
%
%% jn_zeros() computes the zeros of a Bessel function Jn(x).
%
%  Licensing:
%
%    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
%    they give permission to incorporate this routine into a user program
%    provided that the copyright is acknowledged.
%
%  Modified:
%
%    13 June 2025
%
%  Author:
%
%    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
%    This version by John Burkardt.
%
%  Reference:
%
%    Shanjie Zhang, Jianming Jin,
%    Computation of Special Functions,
%    Wiley, 1996,
%    ISBN: 0-471-11963-6,
%    LC: QA351.C45.
%
%  Input:
%
%    integer N: the order of the Bessel functions.
%
%    integer NT: the number of zeros.
%
%  Output:
%
%    real RJ0(NT): the first NT zeros of Jn(x).
%
  rj0 = zeros ( nt, 1 );

  if ( n <= 20 )
    x = 2.82141 + 1.15859 * n;
  else
    x = n + 1.85576 * n^0.33333 + 1.03315 / n^0.33333;
  end

  l = 0;

  while ( true )

    x0 = x;
    [ bjn, djn ] = jyndd ( n, x );
    x = x - bjn / djn;

    if ( 1.0E-09 < abs ( x - x0 ) )
      continue;
    end

    l = l + 1;
    rj0(l) = x;
    x = x + 3.1416 + ( 0.0972 + 0.0679 * n - 0.000354 * n^2 ) / l;

    if ( nt <= l )
      break;
    end

  end

  return
end

