Your technique should have worked, but if you don't know which expansions to do first you can get yourself in a tangle of algebra and make silly mistakes that bring the whole thing crashing down.
The way I reasoned was, well, I have four numbers multiplied together, and I want it to be two numbers of the same size multiplied together. So I'll try multiplying the big one with the small one, and the two middle ones.
$$p(p+1)(p+2)(p+3) + 1 = (p^2 + 3p)(p^2 + 3p + 2) + 1$$
Now those terms are nearly the same. How can we force them together? I'm going to use the basic but sometimes-overlooked fact that $xy = (x+1)y - y$, and likewise $x(y + 1) = xy + x$.
$$\begin{align*}
(p^2 + 3p)(p^2 + 3p + 2) + 1 &= (p^2 + 3p + 1)(p^2 + 3p + 2) - (p^2 + 3p + 2) + 1 \\
&= (p^2 + 3p + 1)(p^2 + 3p + 1) + (p^2 + 3p + 1) - (p^2 + 3p + 2) + 1 \\
&= (p^2 + 3p + 1)^2
\end{align*}$$
Tada.
Suppose there are $n$ and $k$ such that $$ (a_1!-1)(a_2!-1)...(a_n!-1)=k^2+16=k^2+4^2\quad (*).$$
Notice that $4$ can not divide the LHS, and consequently $\mathrm{gcd}(k,4)=1$. Then by step $4$ here, every factor of $k^2+4^2$ is a sum of two squares. That is, for all $i=1,\ldots,n$, there are integers $a$ and $b$ such that $$a_i!-1=a^2+b^2.$$
If $a_i\geq 4$, then $a_i!-1\equiv -1\pmod 4$, but $a^2+b^2\not\equiv -1 \pmod 4$. Hence we must have $a_i\leq 3$ for all $i=1,\ldots,n$.
It follows that the LHS in $(*)$ is either zero, which is impossible, or equal to $$(3!-1)^m=5^m$$
for some $0\leq m \leq n$.
Now let's solve $$ k^2+4^2=5^m.$$
- If $m$ is even, say $m=2c$, then $k^2+4^2=5^{2c}$ yields $2^4=4^2=(5^c-k)(5^c+k)$. Consequently $5^c-k=2^{u}$ and $5^c+k=2^v$ for some $u,v\geq 0$ with $u+v=4$. This is easy to solve, by just discussing cases.
- If $m$ is odd, say $m=2c+1$, then $k^2+4^2=5\cdot 5^{2c}$ yields $$k^2=((2\omega -1) 5^c-4)((2\omega -1) 5^c+4),$$
where $w=\frac{1+\sqrt{5}}{2}$.
The ring $\mathbb{Z}[\omega]$, which is the ring of integers of $\mathbb{Q}(\sqrt{5})$, is a UFD; see oeis.org/A003172. Moreover $(2\omega -1) 5^c-4$ and $(2\omega -1) 5^c+4$ are coprime (because if $p$ is irreducible and divides $(2\omega -1) 5^c-4$ and $(2\omega -1) 5^c+4$, then $p\bar{p}\mid 8$ and so $p\bar{p}=2\mid k^2$, which is impossible). Hence both factors are squares in $\mathbb{Z}[\omega]$, that is, $$(2\omega -1) 5^c-4=(e+\omega f)^2$$
for some $e,f\in \mathbb{Z}$. But using the fact that $\omega^2= 1+\omega$, we can see that the last equation has no solutions in $e,f\in \mathbb{Z}$.
Consequently all the solutions, if any (as I didn't do the computation), come from the case where $m$ is even.
Best Answer
Strategy: try to find two perfect squares such that the expression $f(n)=n^4+6n^3+11n^2+3n+31$ is strictly between them for all but finitely many $n$.
$(n^2+3n+1)^2$ is important to compute, as you also pointed out. It is $n^4+6n^3+11n^2+6n+1$.
So our impression is that $f(n)$ is close to this, and that if $f(n)$ is a perfect square, it should be the square of $n^2+3n+1$.
We prove this for all but finitely many $n$. More precisely, we prove that for all but finitely many $n$, we have $(n^2+3n)^2< f(n)< (n^2+3n+2)^2$. First part: $(n^2+3n)^2= n^4+6n^3+9n^2< n^4+6n^3+11n^2+3n+31 \Leftrightarrow 0< 2n^2+3n+31$. This holds for all $n\in \mathbb{Z}$ (the discriminant is negative). Second part: $n^4+6n^3+11n^2+3n+31< (n^2+3n+2)^2 = n^4 + 6n^3 +13n^2 + 12n + 4 \Leftrightarrow$ $0< 2n^2+9n-27$.
This inequality fails iff $-6\leq n\leq 2$. So unfortunately, you have to check if $f(n)$ is a perfect square for these nine values. (Easy calculation.) Once you are done with that, if $n$ is not among these nine values, then $f(n)$ is strictly between $(n^2+3n)^2$ and $(n^2+3n+2)^2$, so it has to be $(n^2+3n+1)^2$. This leads to a unary equation that you have already solved.