diff --git a/fvn_fnlib/Makefile b/fvn_fnlib/Makefile
index c9def1c..988177a 100644
--- a/fvn_fnlib/Makefile
+++ b/fvn_fnlib/Makefile
@@ -61,7 +61,9 @@ zasin.o zatan2.o zatanh.o zatan.o  \
 zbeta.o zcbrt.o zcosh.o zcot.o  \
 zexprl.o zgamma.o zgamr.o zlbeta.o  \
 zlngam.o zlnrel.o zlog10.o zpsi.o  \
-zsinh.o ztanh.o ztan.o
+zsinh.o ztanh.o ztan.o besyn.o \
+besjn.o dbesyn.o dbesjn.o beskn.o \
+besin.o dbeskn.o dbesin.o
 
 lib:$(library)
 
diff --git a/fvn_fnlib/besin.f90 b/fvn_fnlib/besin.f90
new file mode 100644
index 0000000..5b6db4f
--- /dev/null
+++ b/fvn_fnlib/besin.f90
@@ -0,0 +1,79 @@
+real(4) function besin(n,x,factor,big)
+    implicit none
+    ! This function compute the rank n Bessel J function
+    ! using recurrence relation :
+    ! In+1(x)=-2n/x * In(x) + In-1(x)
+    !
+    ! Two optional parameters :
+    ! factor : an integer that is used in Miller's algorithm to determine the
+    ! starting point of iteration. Default value is 40, an increase of this value
+    ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
+    ! big : a real that determine the threshold for taking anti overflow counter measure
+    ! default value is 1e10
+    ! 
+    integer :: n
+    real(4) :: x
+    integer, optional :: factor
+    real(4), optional :: big
+
+    integer :: tfactor
+    real(4) :: tbig,tsmall
+    real(4) :: two_on_x,binm1,bin,binp1,absx
+    integer :: i,start
+    real(4), external :: besi0,besi1
+
+    ! Initialization of optional parameters
+    tfactor=40
+    if(present(factor)) tfactor=factor
+    tbig=1e10
+    if(present(big)) tbig=big
+    tsmall=1./tbig
+
+    if (n==0) then
+        besin=besi0(x)
+        return
+    end if
+    if (n==1) then
+        besin=besi1(x)
+        return
+    end if
+    if (n < 0) then
+        write(*,*) "Error in besin, n must be >= 0"
+        stop
+    end if
+
+    absx=abs(x)
+    if (absx == 0.) then
+        besin=0.
+    else
+        ! We use Miller's Algorithm
+        ! as upward reccurence is unstable.
+        ! This is adapted from Numerical Recipes
+        ! Principle : use of downward recurrence from an arbitrary
+        ! higher than n value with an arbitrary seed,
+        ! and then use the normalization formula :
+        ! 1=I0-2I2+2I4-2I6+.... however it is easier to use a
+        ! call to besi0
+        two_on_x=2./absx
+        start=2*((n+int(sqrt(float(n*tfactor))))/2) ! even start
+        binp1=0.
+        bin=1.
+        do i=start,1,-1
+            ! begin downward rec
+            binm1=two_on_x*bin*i+binp1
+            binp1=bin
+            bin=binm1
+            ! Action to prevent overflow
+            if (abs(bin) > tbig) then
+                bin=bin*tsmall
+                binp1=binp1*tsmall
+                besin=besin*tsmall
+            end if
+            if (i==n) besin=binp1
+        end do
+        besin=besin*besi0(x)/bin
+    end if
+    ! if n is odd and x <0
+    if ((x<0.) .and. (mod(n,2)==1)) besin=-besin
+
+end function
diff --git a/fvn_fnlib/besjn.f90 b/fvn_fnlib/besjn.f90
new file mode 100644
index 0000000..b4ebb60
--- /dev/null
+++ b/fvn_fnlib/besjn.f90
@@ -0,0 +1,98 @@
+real(4) function besjn(n,x,factor,big)
+    implicit none
+    ! This function compute the rank n Bessel J function
+    ! using recurrence relation :
+    ! Jn+1(x)=2n/x * Jn(x) - Jn-1(x)
+    !
+    ! Two optional parameters :
+    ! factor : an integer that is used in Miller's algorithm to determine the
+    ! starting point of iteration. Default value is 40, an increase of this value
+    ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
+    ! big : a real that determine the threshold for taking anti overflow counter measure
+    ! default value is 1e10
+    ! 
+    integer :: n
+    real(4) :: x
+    integer, optional :: factor
+    real(4), optional :: big
+
+    integer :: tfactor
+    real(4) :: tbig,tsmall,som
+    real(4),external :: besj0,besj1
+    real(4) :: two_on_x,bjnm1,bjn,bjnp1,absx
+    integer :: i,start
+    logical :: iseven
+
+    ! Initialization of optional parameters
+    tfactor=40
+    if(present(factor)) tfactor=factor
+    tbig=1e10
+    if(present(big)) tbig=big
+    tsmall=1./tbig
+
+    if (n==0) then
+        besjn=besj0(x)
+        return
+    end if
+    if (n==1) then
+        besjn=besj1(x)
+        return
+    end if
+    if (n < 0) then
+        write(*,*) "Error in besjn, n must be >= 0"
+        stop
+    end if
+
+    absx=abs(x)
+    if (absx == 0.) then
+        besjn=0.
+    else if (absx > float(n)) then
+    ! For x > n upward reccurence is stable
+        two_on_x=2./absx
+        bjnm1=besj0(absx)
+        bjn=besj1(absx)
+        do i=1,n-1
+            bjnp1=two_on_x*bjn*i-bjnm1
+            bjnm1=bjn
+            bjn=bjnp1
+        end do
+        besjn=bjnp1
+    else
+        ! For x <= n we use Miller's Algorithm
+        ! as upward reccurence is unstable.
+        ! This is adapted from Numerical Recipes
+        ! Principle : use of downward recurrence from an arbitrary
+        ! higher than n value with an arbitrary seed,
+        ! and then use the normalization formula :
+        ! 1=J0+2J2+2J4+2J6+....
+        two_on_x=2./absx
+        start=2*((n+int(sqrt(float(n*tfactor))))/2) ! even start
+        som=0.
+        iseven=.false.
+        bjnp1=0.
+        bjn=1.
+        do i=start,1,-1
+            ! begin downward rec
+            bjnm1=two_on_x*bjn*i-bjnp1
+            bjnp1=bjn
+            bjn=bjnm1
+            ! Action to prevent overflow
+            if (abs(bjn) > tbig) then
+                bjn=bjn*tsmall
+                bjnp1=bjnp1*tsmall
+                besjn=besjn*tsmall
+                som=som*tsmall
+            end if
+            if (iseven) then
+                som=som+bjn
+            end if
+            iseven= .not. iseven
+            if (i==n) besjn=bjnp1
+        end do
+        som=2.*som-bjn
+        besjn=besjn/som
+    end if
+    ! if n is odd and x <0
+    if ((x<0.) .and. (mod(n,2)==1)) besjn=-besjn
+
+end function
diff --git a/fvn_fnlib/beskn.f90 b/fvn_fnlib/beskn.f90
new file mode 100644
index 0000000..cf5f906
--- /dev/null
+++ b/fvn_fnlib/beskn.f90
@@ -0,0 +1,41 @@
+real(4) function beskn(n,x)
+    implicit none
+    ! This function compute the rank n Bessel Y function
+    ! using recurrence relation :
+    ! Kn+1(x)=2n/x * Kn(x) + Kn-1(x)
+    !
+    integer :: n
+    real(4) :: x
+
+    real(4),external :: besk0,besk1
+    real(4) :: two_on_x,bknm1,bkn,bktmp
+    integer :: i
+
+    if (n==0) then
+        beskn=besk0(x)
+        return
+    end if
+    if (n==1) then
+        beskn=besk1(x)
+        return
+    end if
+
+    if (n < 0) then
+        write(*,*) "Error in beskn, n must be >= 0"
+        stop
+    end if
+    if (x <= 0.) then
+        write(*,*) "Error in beskn, x must be strictly positive"
+    end if
+
+    two_on_x=2./x
+    bknm1=besk0(x)
+    bkn=besk1(x)
+
+    do i=1,n-1
+        bktmp=two_on_x*bkn*i+bknm1
+        bknm1=bkn
+        bkn=bktmp
+    end do
+    beskn=bktmp
+end function
diff --git a/fvn_fnlib/besyn.f90 b/fvn_fnlib/besyn.f90
new file mode 100644
index 0000000..9188969
--- /dev/null
+++ b/fvn_fnlib/besyn.f90
@@ -0,0 +1,41 @@
+real(4) function besyn(n,x)
+    implicit none
+    ! This function compute the rank n Bessel Y function
+    ! using recurrence relation :
+    ! Yn+1(x)=2n/x * Yn(x) - Yn-1(x)
+    !
+    integer :: n
+    real(4) :: x
+
+    real(4),external :: besy0,besy1
+    real(4) :: two_on_x,bynm1,byn,bytmp
+    integer :: i
+
+    if (n==0) then
+        besyn=besy0(x)
+        return
+    end if
+    if (n==1) then
+        besyn=besy1(x)
+        return
+    end if
+
+    if (n < 0) then
+        write(*,*) "Error in besyn, n must be >= 0"
+        stop
+    end if
+    if (x <= 0.) then
+        write(*,*) "Error in besyn, x must be strictly positive"
+    end if
+
+    two_on_x=2./x
+    bynm1=besy0(x)
+    byn=besy1(x)
+
+    do i=1,n-1
+        bytmp=two_on_x*byn*i-bynm1
+        bynm1=byn
+        byn=bytmp
+    end do
+    besyn=bytmp
+end function
diff --git a/fvn_fnlib/dbesin.f90 b/fvn_fnlib/dbesin.f90
new file mode 100644
index 0000000..a71321f
--- /dev/null
+++ b/fvn_fnlib/dbesin.f90
@@ -0,0 +1,80 @@
+real(8) function dbesin(n,x,factor,big)
+    implicit none
+    ! This function compute the rank n Bessel J function
+    ! using recurrence relation :
+    ! In+1(x)=-2n/x * In(x) + In-1(x)
+    !
+    ! Two optional parameters :
+    ! factor : an integer that is used in Miller's algorithm to determine the
+    ! starting point of iteration. Default value is 40, an increase of this value
+    ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
+    ! big : a real that determine the threshold for taking anti overflow counter measure
+    ! default value is 1e10
+    ! 
+    integer :: n
+    real(8) :: x
+    integer, optional :: factor
+    real(8), optional :: big
+
+    integer :: tfactor
+    real(8) :: tbig,tsmall
+    real(8) :: two_on_x,binm1,bin,binp1,absx
+    integer :: i,start
+    real(8), external :: dbesi0,dbesi1
+
+    ! Initialization of optional parameters
+    tfactor=40
+    if(present(factor)) tfactor=factor
+    tbig=1e10
+    if(present(big)) tbig=big
+    tsmall=1./tbig
+
+    if (n==0) then
+        dbesin=dbesi0(x)
+        return
+    end if
+    if (n==1) then
+        dbesin=dbesi1(x)
+        return
+    end if
+
+    if (n < 0) then
+        write(*,*) "Error in dbesin, n must be >= 0"
+        stop
+    end if
+
+    absx=abs(x)
+    if (absx == 0.) then
+        dbesin=0.
+    else
+        ! We use Miller's Algorithm
+        ! as upward reccurence is unstable.
+        ! This is adapted from Numerical Recipes
+        ! Principle : use of downward recurrence from an arbitrary
+        ! higher than n value with an arbitrary seed,
+        ! and then use the normalization formula :
+        ! 1=I0-2I2+2I4-2I6+.... however it is easier to use a
+        ! call to besi0
+        two_on_x=2./absx
+        start=2*((n+int(sqrt(float(n*tfactor))))/2) ! even start
+        binp1=0.
+        bin=1.
+        do i=start,1,-1
+            ! begin downward rec
+            binm1=two_on_x*bin*i+binp1
+            binp1=bin
+            bin=binm1
+            ! Action to prevent overflow
+            if (abs(bin) > tbig) then
+                bin=bin*tsmall
+                binp1=binp1*tsmall
+                dbesin=dbesin*tsmall
+            end if
+            if (i==n) dbesin=binp1
+        end do
+        dbesin=dbesin*dbesi0(x)/bin
+    end if
+    ! if n is odd and x <0
+    if ((x<0.) .and. (mod(n,2)==1)) dbesin=-dbesin
+
+end function
diff --git a/fvn_fnlib/dbesjn.f90 b/fvn_fnlib/dbesjn.f90
new file mode 100644
index 0000000..da26cd1
--- /dev/null
+++ b/fvn_fnlib/dbesjn.f90
@@ -0,0 +1,98 @@
+real(8) function dbesjn(n,x,factor,big)
+    implicit none
+    ! This function compute the rank n Bessel J function
+    ! using recurrence relation :
+    ! Jn+1(x)=2n/x * Jn(x) - Jn-1(x)
+    !
+    ! Two optional parameters :
+    ! factor : an integer that is used in Miller's algorithm to determine the
+    ! starting point of iteration. Default value is 40, an increase of this value
+    ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
+    ! big : a real that determine the threshold for taking anti overflow counter measure
+    ! default value is 1e10
+    ! 
+    integer :: n
+    real(8) :: x
+    integer, optional :: factor
+    real(8), optional :: big
+
+    integer :: tfactor
+    real(8) :: tbig,tsmall,som
+    real(8),external :: dbesj0,dbesj1
+    real(8) :: two_on_x,bjnm1,bjn,bjnp1,absx
+    integer :: i,start
+    logical :: iseven
+
+    ! Initialization of optional parameters
+    tfactor=40
+    if(present(factor)) tfactor=factor
+    tbig=1d10
+    if(present(big)) tbig=big
+    tsmall=1./tbig
+
+    if (n==0) then
+        dbesjn=dbesj0(x)
+        return
+    end if
+    if (n==1) then
+        dbesjn=dbesj1(x)
+        return
+    end if
+    if (n < 0) then
+        write(*,*) "Error in dbesjn, n must be >= 0"
+        stop
+    end if
+
+    absx=abs(x)
+    if (absx == 0.) then
+        dbesjn=0.
+    else if (absx > float(n)) then
+    ! For x > n upward reccurence is stable
+        two_on_x=2./absx
+        bjnm1=dbesj0(absx)
+        bjn=dbesj1(absx)
+        do i=1,n-1
+            bjnp1=two_on_x*bjn*i-bjnm1
+            bjnm1=bjn
+            bjn=bjnp1
+        end do
+        dbesjn=bjnp1
+    else
+        ! For x <= n we use Miller's Algorithm
+        ! as upward reccurence is unstable.
+        ! This is adapted from Numerical Recipes
+        ! Principle : use of downward recurrence from an arbitrary
+        ! higher than n value with an arbitrary seed,
+        ! and then use the normalization formula :
+        ! 1=J0+2J2+2J4+2J6+....
+        two_on_x=2./absx
+        start=2*((n+int(sqrt(float(n*tfactor))))/2) ! even start
+        som=0.
+        iseven=.false.
+        bjnp1=0.
+        bjn=1.
+        do i=start,1,-1
+            ! begin downward rec
+            bjnm1=two_on_x*bjn*i-bjnp1
+            bjnp1=bjn
+            bjn=bjnm1
+            ! Action to prevent overflow
+            if (abs(bjn) > tbig) then
+                bjn=bjn*tsmall
+                bjnp1=bjnp1*tsmall
+                dbesjn=dbesjn*tsmall
+                som=som*tsmall
+            end if
+            if (iseven) then
+                som=som+bjn
+            end if
+            iseven= .not. iseven
+            if (i==n) dbesjn=bjnp1
+        end do
+        som=2.*som-bjn
+        dbesjn=dbesjn/som
+    end if
+    ! if n is odd and x <0
+    if ((x<0.) .and. (mod(n,2)==1)) dbesjn=-dbesjn
+
+end function
diff --git a/fvn_fnlib/dbeskn.f90 b/fvn_fnlib/dbeskn.f90
new file mode 100644
index 0000000..083c428
--- /dev/null
+++ b/fvn_fnlib/dbeskn.f90
@@ -0,0 +1,41 @@
+real(8) function dbeskn(n,x)
+    implicit none
+    ! This function compute the rank n Bessel Y function
+    ! using recurrence relation :
+    ! Kn+1(x)=2n/x * Kn(x) + Kn-1(x)
+    !
+    integer :: n
+    real(8) :: x
+
+    real(8),external :: dbesk0,dbesk1
+    real(8) :: two_on_x,bknm1,bkn,bktmp
+    integer :: i
+
+    if (n==0) then
+        dbeskn=dbesk0(x)
+        return
+    end if
+    if (n==1) then
+        dbeskn=dbesk1(x)
+        return
+    end if
+
+    if (n < 0) then
+        write(*,*) "Error in dbeskn, n must be >= 0"
+        stop
+    end if
+    if (x <= 0.) then
+        write(*,*) "Error in dbeskn, x must be strictly positive"
+    end if
+
+    two_on_x=2./x
+    bknm1=dbesk0(x)
+    bkn=dbesk1(x)
+
+    do i=1,n-1
+        bktmp=two_on_x*bkn*i+bknm1
+        bknm1=bkn
+        bkn=bktmp
+    end do
+    dbeskn=bktmp
+end function
diff --git a/fvn_fnlib/dbesyn.f90 b/fvn_fnlib/dbesyn.f90
new file mode 100644
index 0000000..41b133f
--- /dev/null
+++ b/fvn_fnlib/dbesyn.f90
@@ -0,0 +1,40 @@
+real(8) function dbesyn(n,x)
+    implicit none
+    ! This function compute the rank n Bessel Y function
+    ! using recurrence relation :
+    ! Yn+1(x)=2n/x * Yn(x) - Yn-1(x)
+    !
+    integer :: n
+    real(8) :: x
+
+    real(8),external :: dbesy0,dbesy1
+    real(8) :: two_on_x,bynm1,byn,bytmp
+    integer :: i
+
+    if (n==0) then
+        dbesyn=dbesy0(x)
+        return
+    end if
+    if (n==1) then
+        dbesyn=dbesy1(x)
+        return
+    end if
+    if (n < 0) then
+        write(*,*) "Error in dbesyn, n must be >= 0"
+        stop
+    end if
+    if (x <= 0.) then
+        write(*,*) "Error in dbesyn, x must be strictly positive"
+    end if
+
+    two_on_x=2./x
+    bynm1=dbesy0(x)
+    byn=dbesy1(x)
+
+    do i=1,n-1
+        bytmp=two_on_x*byn*i-bynm1
+        bynm1=byn
+        byn=bytmp
+    end do
+    dbesyn=bytmp
+end function
diff --git a/fvn_fnlib/fvn_fnlib.f90 b/fvn_fnlib/fvn_fnlib.f90
index bb6dee3..20079a1 100644
--- a/fvn_fnlib/fvn_fnlib.f90
+++ b/fvn_fnlib/fvn_fnlib.f90
@@ -733,6 +733,62 @@ interface bsk1e
       end function dbsk1e
 end interface bsk1e
 
+! nth order J
+interface bsjn
+      real(4) function besjn(n,x,factor,big)
+            integer(4) :: n
+            real(4) :: x
+            integer(4), optional :: factor
+            real(4), optional :: big
+      end function besjn
+      real(8) function dbesjn(n,x,factor,big)
+            integer(4) :: n
+            real(8) :: x
+            integer(4), optional :: factor
+            real(8), optional :: big
+      end function dbesjn
+end interface bsjn
+
+! nth order Y
+interface bsyn
+      real(4) function besyn(n,x)
+            integer(4) :: n
+            real(4) :: x
+      end function besyn
+      real(8) function dbesyn(n,x)
+            integer(4) :: n
+            real(8) :: x
+      end function dbesyn
+end interface bsyn
+
+! nth order I
+interface bsin
+      real(4) function besin(n,x,factor,big)
+            integer(4) :: n
+            real(4) :: x
+            integer(4), optional :: factor
+            real(4), optional :: big
+      end function besin
+      real(8) function dbesin(n,x,factor,big)
+            integer(4) :: n
+            real(8) :: x
+            integer(4), optional :: factor
+            real(8), optional :: big
+      end function dbesin
+end interface bsin
+
+! nth order K
+interface bskn
+      real(4) function beskn(n,x)
+            integer(4) :: n
+            real(4) :: x
+      end function beskn
+      real(8) function dbeskn(n,x)
+            integer(4) :: n
+            real(8) :: x
+      end function dbeskn
+end interface bskn
+
 !!!!!!!!!!!!!!!!!!!!!
 ! MISSING BSJNS
 ! MISSING BSINS