program test_bsjn
  use fvn_fnlib
  implicit none
  integer :: n,npoints,i

  real(kind=sp_kind) :: xmin,xmax,xstep,x,y
  real(kind=sp_kind),dimension(200,0:5) :: bes

  open(2,file='bsjn.dat')
  
  xmin=-20.
  xmax=20.
  npoints=200
  xstep=(xmax-xmin)/dble(npoints)
  do i=1,npoints
        x=xmin+i*xstep
        do n=0,5
          bes(i,n)=bsjn(n,x)
        end do
        write(2,'(7e22.14)') x,bes(i,0:5)
  end do
  close(2)
  
  end program