I think you have mistakes. Please, see the condition \[ a_{n+1}=a_n + ( a_1+1 ) ( 2 a_2 + 1) ( 3a_3+1)\cdots (n - 1)a_{n - 1}+1)((n^2 + n)a_n+2n+1)) .\]

Yeah, there is a parenthesis before the $a_{n-1}$'s coefficient. I've just edited it.

something is still wrong with the last term...

Yeah, although it looks like wrong, this is the problem written on the paper.

Does anybody has a solution to this strange problem?

Any body solution?

It's enough to prove that $b_n=na_n+5n^2+1$ is always a perfect square. For $n=1,2,3,4,5$ it equals $1^2,3^2,20^2,501^2,313026^2$, respectively. So i think it's true

It surely does!!! But I don't know how to solve it.

Let $b_n=na_n+1$ and $B_n=\prod_{i=1}^nb_i$. Obviously $2|B_1|B_n$. Also $2|a_n,n>1$ because $a_{n+1}=a_n\mod2,n>1$, The equation becomes $a_{n+1}=a_n+B_{n-1}((n+1)b_n+n)=a_n+(n+1)B_n+nB_{n-1}$. Multiply by $n+1$ and plus $1$, $b_{n+1}=(n+1)a_n+1+(n+1)^2B_n+(n^2+n)B_{n-1}$ Multiply by $B_n$, $B_{n+1}=(n+1)^2B_n^2+(n^2+n)B_nB_{n-1}+(n+1)a_nB_n+B_n \\=[(n+1)B_n+\frac n2B_{n-1}+\frac{a_n}2]^2+B_n-\frac14(nB_{n-1}+a_n)^2$. Hence, $B_{n+1}-\frac14((n+1)B_n+a_{n+1})^2 \\=B_{n+1}-\frac14(2(n+1)B_n+a_n+nB_{n-1})^2 \\=B_{n+1}-[(n+1)B_n+\frac n2B_{n-1}+\frac{a_n}2]^2 \\=B_n-\frac14(nB_{n-1}+a_n)^2$ Now $B_n-\frac14(nB_{n-1}+a_n)^2=\cdots=B_2-\frac14(2B_1+a_2)^2=-5,\forall n>1$ In other words, we can express $B_n$ as $x^2-5$. Since $na_n+1=b_n|B_n$, we are done!

WOW it is FANTASTIC!!