For Fermi-Dirac Integrals
$$\mathcal{F}_{j}(x)=\frac{1}{\Gamma(j+1)}\int_{0}^{\infty}\frac{t^{j}}{e^{t-x}+1}~\mathrm{d}t,$$
we know that for $x\gg1,$ we can approximate
$$\mathcal{F}_{0}(x) \approx x$$
and
$$\mathcal{F}_{1}(x) \approx\frac{1}{2}x^{2}.$$
However, using numerical integration for obtaining the integrals, the function $2\mathcal{F}_{0}(x)-\mathcal{F}_{0}(x)^{2}$ converges to a constant threshold.
How can that be?