zeroj(m_order,n_zero)

function z = zeroj(m_order,n_zero)
% Function which returns the zeros of the Bessel function J(x)
%   z = zeroj(m_order,n_zeros)
% Inputs
%   m_order - Order of the Bessel function
%   n_zero  - Number of the zero
% Output
%   z - The "n_zero th" zero of the Bessel function
%% Use asymtotic formula for initial guess
beta = (n_zero + 0.5*m_order - 0.25)*pi;
mu = 4*m_order^2;
z = beta - (mu-1)/(8*beta) - 4*(mu-1)*(7*mu-31)/(3*(8*beta)^3);
%fprintf('Initial guess is %15.10g \n',z);
for i=1:5
  jj = bess(m_order+1,z);          % Use the recursion relation
  deriv = -jj(m_order+2) + ...     % to evaluate derivative
         m_order/z * jj(m_order+1);
  z = z - jj(m_order+1)/deriv;     % Newton's root finding  
%  fprintf('Iteration %g; value = %15.10g \n',i,z);
end
return;