Help with QAGI from QUADPACK



While testing QAGI from QUADPACK I ran across a problem that seems to
expose a problem in QAGI. I'm wondering if anyone here has seen
similar behavior and/or knows of a fix to QAGI.
The code below attemps to integrate the Normal PDF from negative
infinity to a bound. When I set the bound to something like 33 (or
less) I get reasonable results. When the bound is set to 34 or larger
the result approaches zero, without any error condition being flagged
by QAGI.
Any help is appreciated. BTW, I really want to stick with QUADPACK, so
suggesting an alternative package is not really an option for me. I'm
more interested in understanding what is going on in QAGI.

I downloaded QAGI from netlib, and I get the same results with both the
Sun Compiler on Solaris and Compaq Visual Fortran on Windows.

Thank you,
--Delta

program test
implicit none
integer :: limit, lenw
parameter (limit=500, lenw=4*limit)
integer :: i,j, neval, ier, last, iwork(limit)
real :: err,y, bound, work(lenw)
real, external :: FCN

bound = 34.0
call qagi(fcn, bound, -1, 0.0001, 0.0, y, err, neval, ier, &
limit, lenw, last, iwork, work)

write(*,*) "QAGI: ans = ", y
write(*,*) "QAGI: ier = ",ier
end program

! Normal PDF
real function FCN(x)
implicit none
real :: rsigsq, scale, sigma, mu, pi
real, intent(in) :: x
intrinsic atan, exp, sqrt
pi=4*atan(1.)
sigma = 1.0
mu = 0.0

rsigsq=1.e0/(sqrt(2.e0)*sigma)
scale=1.e0/sqrt(2.e0*pi)
FCN=(scale/sigma)*exp(-((x-mu)*rsigsq)**2 )
! write(*,*) x, FCN
RETURN
end function

.