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