Consider Fermi-Dirac-like integral, with $x$ positive real number and $a$ positive real parameter
$$I(a) = \int\frac{x^3\cdot\exp{\left(\dfrac{x-1}{a}\right)}}{\left[1+\exp{\left(\dfrac{x-1}{a}\right)}\right]^2}~\mathrm{d}x$$.
Mathematica spits out ugly solution with polylog.
In extreme limit $a \to 0$, integrand is delta function and $I(a)$ is step function, allowing simple solution.
How to get approximate solution for $I(a)$ in regime $a \ll 1$ (and not brutal $a \to 0$ making use of delta function property)?
We have the following approximation:
$$ \int_{-\infty}^{x} t^3 \frac{e^{-(t-1)/a}}{(1 + e^{-(t-1)/a})^2} \, \mathrm{d}t = \frac{a}{1 + e^{-(x-1)/a}} + \mathcal{O}(a^2),$$
where the implicit constant for the asymptotic notation is uniform in $x$ and $a$. As a result,
$$ \lim_{a\to 0^+} \frac{1}{a}\int_{-\infty}^{x} t^3 \frac{e^{-(t-1)/a}}{(1 + e^{-(t-1)/a})^2} \, \mathrm{d}t = \begin{cases} 1, & x > 1 \\ 1/2, & x = 1 \\ 0, & x < 1 \end{cases}. $$
The idea is extremely simple. Write $u = (t-1)/2a$ and notice that
\begin{align*} \int_{-\infty}^{x} t^3 \frac{e^{-(t-1)/a}}{(1 + e^{-(t-1)/a})^2} \, \mathrm{d}t &= \int_{-\infty}^{\frac{x-1}{2a}} a (2au + 1)^3 \frac{2 e^{-2u}}{(1 + e^{-2u})^2} \, \mathrm{d}u \\ &= \int_{-\infty}^{\frac{x-1}{2a}} (a + 6a^2u + \cdots) \frac{2 e^{-2u}}{(1 + e^{-2u})^2} \, \mathrm{d}u. \tag{*} \end{align*}
Using a simple estimate $ \frac{e^{-2u}}{(1 + e^{-2u})^2} \leq e^{-2|u|}$, it follows that
$$ \int_{-\infty}^{\infty} |u|^s \frac{e^{-2u}}{(1 + e^{-2u})^2} \, \mathrm{d}u < \infty \qquad \forall s \geq 0.$$
This allows us to safely ignore the higher order terms in $\text{(*)}$, as they only contribute to $\mathcal{O}(a^2)$ term. From this, we can write
$$ \text{(*)} = a \int_{-\infty}^{\frac{x-1}{2a}} \frac{2 e^{-2u}}{(1 + e^{-2u})^2} \, \mathrm{d}u + \mathcal{O}(a^2). $$
Evaluating the last integral gives the desired result. We can also identify higher order terms, but they involve polylogarithms.