program test_bskn
use fvn_fnlib
use Kind_Definition
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='bskn.dat')

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

! with gnuplot -> set yrange [-1:1]

end program