MS,
Looking at your code finds two problems:
1) Your calculation is numerically unstable. You can see this by
looking at taum(1.0,...) in double precision and in extended-precision,
bfloat(taum(1.0,...)),fpprec:30. If you look at the analytical form,
taum(x,...) you will see why. The denominator is of the form
log(x+1)-x/(x+1). Since your x's are very near zero (1.0e-11), you lose
most precision in calculating log(x+1). You can solve this problem
either by rewriting the denominator as x^2/2 (verify using Maxima's
Taylor command if you like) or by calculating everything in extended
precision.
2) plot2d has a bug in scaling functions with large ranges. Bug report
https://sourceforge.net/tracker/index.php?func=detail&aid=694770&group_id=4933&atid=104933.
You will have to scale "by hand", by multiplying by
1.0e-11 or whatever.
Note by the way that for parameters in the range you're looking at, taum
is nearly exactly k*energy^(-1/2), for a value of 'k' which I'll leave
you the pleasure of rederiving.
-s