It's an integral inequality I find delightful let me share it :
$$I=\int_{0}^{3}\frac{x!}{x^{2}+1}dx>\frac{1+\sqrt{5}}{2}$$
To show it I have managed to use Stirling's approximation as we have :
$$\int_{0}^{3}\frac{\sqrt{2\pi x}e^{\frac{1}{12x+1}}\left(\frac{x}{e}\right)^{x}}{x^{2}+1}dx<I$$
I have also tried Ramanujan inequality but I haven't the skill to evaluate the integral .
Update :
we have the sufficient approximation for $x\in[2,3]$ :
$$\sqrt{2\pi x}\left(1+\left(\frac{x}{e}-1\right)x+\frac{1}{2}\left(\frac{x}{e}-1\right)^{2}x\left(x-1\right)+\frac{1}{6}\left(\frac{x}{e}-1\right)^{3}x\left(x-1\right)\left(x-2\right)\right)e^{\left(\frac{1}{12x+\frac{2}{7}}\right)}-x!<0$$
Where we have used binomial series around $x=1$.
We can also perform a little conjecture for $x\in[0,1]$ :
$$1+x^{\frac{9}{10}}\left(\frac{\pi-e+\frac{1}{2}}{2}x^{\frac{92}{100}}-\frac{\pi-e+\frac{1}{2}}{2}\right)+\frac{1}{32}x^{2}\left(1-x\right)^{2}+\frac{1}{6}x^{4}\left(1-x\right)^{4}\leq^? x!$$
Update 2 :
We have :
$$\int_{0}^{1}\ln\left(x!+1\right)dx+1<\int_{0}^{1}\ln\left(\frac{x!}{x^{2}+1}\right)dx+2$$
$$\int_{0}^{1}\frac{\ln^{2}\left(x!+1\right)}{1+\frac{1}{2}+\sqrt{\frac{5}{4}}}dx<\int_{0}^{1}\ln^{2}\left(\frac{x!}{x^{2}+1}\right)dx$$
$$1-\frac{\sqrt{\frac{5}{4}}-\frac{1}{2}}{2}\int_{0}^{1}\left(\ln\left(x!+1\right)\right)^{3}dx<\int_{0}^{1}\left(\ln\left(\frac{x!}{x^{2}+1}\right)\right)^{3}dx+1$$
$$\left(\int_{0}^{1}\left(\ln\left(x!+1\right)\right)^{4}dx\right)\left(\sqrt{\frac{5}{4}}-\frac{1}{2}\right)^{2}\cdot\frac{20}{29}<\int_{0}^{1}\left(\ln\left(\frac{x!}{x^{2}+1}\right)\right)^{4}dx$$
Now we can use :
$$x=1+\ln(x)+\cdots+\frac{\left(\ln x\right)^{n}}{n!}+\cdots$$
We have also $x\in[0,1]$:
$$x!\simeq \cos\left(\frac{9}{10}\sqrt{\frac{1}{4}-\left(-\frac{1}{2}+x\right)^{2}}+\frac{1}{7.25}\sqrt{\left(x\right)\left(1-x\right)^{4}}+\frac{1}{60}\sqrt{x\left(1-\frac{2x}{x+1}\right)}+\frac{1}{10}\sqrt{x}\left(1-x\right)\left(x+1-\sqrt{2}\right)^{2}\right)$$
How to show it analytically ?
Reference :
Evaluation by Wolfram alpha give $I \approx 1.61815$ and $\phi \approx 1.61803$. It makes a diffrence of the order of $10^{-4}$, thus you need very sharp bounds to prove this inequality formally. I think that the most simple (long but simple) would be to use a formal calculus software like Mapple and to use the fact that for any $a < b$ and any convex function (your integrand is convex) on $[a,b]$ and any $c \in [a,b]$, we have, $$ (1) \quad \forall a \leqslant x \leqslant c, f(x) \geqslant \frac{f(b) - f(c)}{b - c}(x - c) + f(c) $$ Let $f : x \mapsto \frac{x!}{x^2 + 1}$. Assume we can compute approximations $f_\leqslant$ and $f_\geqslant$ of $f$ such that $f_\leqslant \leqslant f \leqslant f_\geqslant$. For any $n \geqslant 0$, \begin{align*} I & = \sum_{k = 0}^{3n - 1} \int_{k/n}^{(k + 1)/n} f(x) \, dx\\ & \geqslant \sum_{k = 0}^{3n - 1} \int_{k/n}^{(k + 1)/n} \frac{f\left(\frac{k + 2}{n}\right) - f\left(\frac{k + 1}{n}\right)}{\frac{k + 2}{n} - \frac{k + 1}{n}}\left(x - \frac{k + 1}{n}\right) + f\left(\frac{k + 1}{n}\right) \, dx\\ & \textrm{ using $(1)$ with $a = k/n$, $b = (k + 2)/n$ and $c = (k + 1)/n$}.\\ & = \sum_{k = 0}^{3n - 1} \int_{k/n}^{(k + 1)/n} n\left(f\left(\frac{k + 2}{n}\right) - f\left(\frac{k + 1}{n}\right)\right)\left(x - \frac{k + 1}{n}\right) + f\left(\frac{k + 1}{n}\right) \, dx\\ & = \sum_{k = 0}^{3n - 1} n\left(f\left(\frac{k + 2}{n}\right) - f\left(\frac{k + 1}{n}\right)\right)\left[\frac{1}{2}\left(x - \frac{k + 1}{n}\right)\right]_{k/n}^{(k + 1)/n} + \frac{1}{n}f\left(\frac{k + 1}{n}\right)\\ & = \frac{1}{n}\sum_{k = 0}^{3n - 1} -\frac{1}{2}\left(f\left(\frac{k + 2}{n}\right) - f\left(\frac{k + 1}{n}\right)\right) + f\left(\frac{k + 1}{n}\right)\\ & = \frac{1}{n}\sum_{k = 0}^{3n - 1} -\frac{1}{2}f\left(\frac{k + 2}{n}\right) + \frac{3}{2}f\left(\frac{k + 1}{n}\right)\\ & \geqslant \frac{1}{n}\sum_{k = 0}^{3n - 1} -\frac{1}{2}f_\geqslant\left(\frac{k + 2}{n}\right) + \frac{3}{2}f_\leqslant\left(\frac{k + 1}{n}\right). \end{align*} Now, the idea is to find sequences of $f_{\leqslant,n} \leqslant f \leqslant f_{\geqslant,n}$ such that $f_{\leqslant,n} \rightarrow f$ and $f_{\geqslant,n} \rightarrow f$ uniformly such that each $f_{\leqslant,n}$ and $f_{\geqslant,n}$ sends rational numbers on rational numbers and is computable at such points. It is probably possible but it will ask a lot of calculations with ugly sums/integrals because of the presence of $x!$ in the definition of $f$.
Finally, if you find a large enough integer $N$ such that, $$ \frac{1}{N}\sum_{k = 0}^{3N - 1} -\frac{1}{2}f_{\geqslant,N}\left(\frac{k + 2}{N}\right) + \frac{3}{2}f_{\leqslant,N}\left(\frac{k + 1}{N}\right) > \frac{1 + \sqrt{5}}{2}, $$ that can be checked formally by a computer by verifying that, $$ \left(-1 + \frac{1}{N}\sum_{k = 0}^{3N - 1} -f_{\geqslant,N}\left(\frac{k + 2}{N}\right) + 3f_{\leqslant,N}\left(\frac{k + 1}{N}\right)\right)^2 > 5, $$ then you have your proof.