function [ bjn, djn ] = jyndd ( n, x )

%*****************************************************************************80
%
%% jyndd() evaluates a Bessel function Jn(x) and derivative.
%
%  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.
%
%    real X, the argument.
%
%  Output:
%
%    real BJN, DJN, the values of Jn(x), Jn'(x).
%
  bj = zeros ( n + 2, 1 );

  for nt = 1 : 900
    mt = floor ( 0.5 * log10 ( 6.28 * nt ) ...
      - nt * log10 ( 1.36 * abs ( x ) / nt ) );
    if ( 20 < mt )
      break;
    end
  end

  m = nt;
  bs = 0.0;
  f0 = 0.0;
  f1 = 1.0E-35;
  su = 0.0;

  for k = m : -1 : 0
 
    f = 2.0 * ( k + 1.0 ) * f1 / x - f0;
    if ( k <= n + 1 )
      bj(k+1) = f;
    end
    if ( mod ( k, 2 ) == 0 )
      bs = bs + 2.0 * f;
      if ( k ~= 0 )
        su = su + ( -1.0 )^( k / 2 ) * f / k;
      end
    end
    f0 = f1;
    f1 = f;
  end

  for k = 0 : n + 1
    bj(k+1) = bj(k+1) / ( bs - f );
  end

  bjn = bj(n+1);
  djn = - bj(n+2) + n * bj(n+1) / x;

  return
end
