Hi there,
I'd like to do numerical integration in my Fortran 90 code, so as a test I tried integrating the function y=x^2 from x=10 to x=12 using the Numerical Recipes routine qtrap. The answer I should get is (12^3)/3-(10^3)/3=242.67, but instead the code returns 22.802.
I traced the problem to line 17 in the trapzd code, s=0.5_sp*(b-a)*sum(func( (/ a,b /) )). Here the code should find sum(fun( (/a,b/) )) =10^2+12^2=244, but instead it finds sum(fun( (/a,b/) ))=3.6E-40. I don't understand why it is getting such a tiny number. Can anyone explain what's going on and how to get qtrap to work correctly?
In case it's useful, I'll copy the relevant lines of code below. Thanks in advance for any help!

Thanks,
Jonell

print *, QTRAP(square,10.,12.)

function square(x) result(y)
real :: x, y
y=x**2
end function square