***********************************
*  ベッセル関数のゼロ点を計算する *
***********************************
      subroutine beszero(nmax,mmax,k)
      integer nmax          ! ベッセル関数のゼロ点の最大個数
      integer mmax          ! ベッセル関数の最大次数
      real bessj
      real a, b, c, d, f, g, lim
      real dx
      real k(0:nmax,mmax)

*-- 二分法の解と近似する条件 ---
      lim=1.0d-6

*-- 二分法の二点を決定するための, 刻み幅 ---
*-- ベッセル関数のゼロ点の間隔はおよそ 3 ごとであるので,
*-- 0.5 ずつ刻めば, まあいいか
      dx=0.5d0
*-- 配列の初期化 ---
      do 1 i=0,nmax
         do 2 j=1,mmax
            k(i,j)=0.0d0
  2   continue
  1   continue

*-- 0 次計算 ---
      k(0,1)=0.0d0
      d=k(0,1)
      do 10 i=1,mmax
         if(i.gt.1)then
            d=k(0,i-1)+dx
         end if
         do while (k(0,mmax).eq.0.0d0)
            a=d
            e=bessj(0,a)
            b=a+dx
            f=bessj(0,b)
            d=d+dx
            do while (e*f.lt.0.0d0)
               c=0.5d0*(a+b)
               g=bessj(0,c)
               if(e*g.lt.0.0d0)then
                  b=c
               else
                  a=c
               end if
               if(abs(g).lt.lim)then
                  k(0,i)=c
                  go to 10
               end if
            end do
        end do 
  10  continue
      if(nmax.eq.0)then
        go to 100
      end if
*-- 1 次以上計算 ---
      do 20 n=1,nmax
      do 21 i=1,mmax
         d=k(n-1,i)+dx
         do while (k(n,mmax).eq.0.0d0)
            a=d
            e=bessj(n,a)
            b=a+dx
            f=bessj(n,b)
            d=d+dx
            do while (e*f.lt.0.0d0)
               c=0.5d0*(a+b)
               g=bessj(n,c)
               if(e*g.lt.0.0d0)then
                  b=c
               else
                  a=c
               end if
               if(abs(g).lt.lim)then
                  k(n,i)=c
                  go to 21
               end if
            end do
         end do
  21  continue
  20  continue
  100 return
      end