Debye function
f(x)=3∫01t3/[exp(tx)−1]dt
x (float or 1d array)
nstep (int) – number of points for integration
zero (float) – approximate the 0 in the integral by this (small!) number