I am trying to compute the Bessel function of the second kind (Bessel_y) using the SLATEC's Amos library available on Netlib. Here is the SLATEC code I use. Below I have pasted my test program that calls SLATEC routine CBESY.

PROGRAM BESSELTEST
IMPLICIT NONE
REAL:: FNU
INTEGER, PARAMETER :: N = 2, KODE = 1
COMPLEX,ALLOCATABLE :: CWRK (:), CY (:)
COMPLEX:: Z, ci
INTEGER :: NZ, IERR

ALLOCATE(CWRK(N), CY(N))
ci = cmplx (0.0, 1.0)
FNU = 0.0e0
Z = CMPLX(0.3e0, 0.4e0)
CALL CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)
WRITE(*,*) 'CY: ', CY
WRITE(*,*) 'IERR: ', IERR
STOP
END PROGRAM

And here is the output of the above program:

CY: ( 5.78591091E-39, 5.80327020E-39) ( 0.0000000 , 0.0000000 )
IERR: 4

Ierr = 4 meaning there is some problem with the input itself. To be precise, the IERR = 4 means the following as per the header info in CBESY.f file:

! IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
! TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
! CANCE BY ARGUMENT REDUCTION

Clearly, CABS(Z) (which is 0.50) or FNU + N - 1 (which is 1.0) are not too large but still the routine CBESY throws the error message number 4 as above. The CY array should have following values for the argument given in above code:

CY(1) = -0.4983 + 0.6700i
CY(2) = -1.0149 + 0.9485i

These values are computed using Matlab. I can't figure out what's the problem when I call CBESY from SLATEC library. Any clues? Much thanks for the suggestions/help. PS: if it is of any help, I used gfortran to compile, link and then create the SLATEC library file ( the .a file ) which I keep in the same directory as my test program above. shell command to execute above code:

gfortran -c BesselTest.f95
gfortran -o a *.o libslatec.a
a

GD.