A challenging system of coupled recursive sequences So I was runing through my old school drafts, and I've just come upon this challenging problem that years ago, one of my former math teacher had let for enthousiastic students to try.
Consider the real sequences $\left(a_n\right)_{n\in\mathbb{N}^*}$ and $\left(b_n\right)_{n\in\mathbb{N}^*}$, defined by :
$$\left|\begin{array}{lll}\displaystyle a_1=-1\\\displaystyle b_1=5\end{array}\right.$$
$$\text{and}$$
$$\forall n\geq2,\left|\begin{array}{lll}\displaystyle a_n=na_{n-1}+\frac{n+1}{n^2+1}b_{n-1}\\\displaystyle b_n=-\frac{n^2+2n+2}{n+2}a_{n-1}-\frac{n^4+3n^3+4n^2+2n}{n^3+2n^2+n+2}b_{n-1}\end{array}\right.$$
And the following questions :


*

*(The most difficult) Find closed form expressions for $\left(a_n\right)_{n\in\mathbb{N}^*}$ and $\left(b_n\right)_{n\in\mathbb{N}^*}$

*(A weaker result, thus easier) Find $\lim\limits_{n\to\infty}\frac{b_n}{n^2a_n}$.


I should say it upfront, there $are$ closed form expressions of $\left(a_n\right)_{n\in\mathbb{N}^*}$,  $\left(b_n\right)_{n\in\mathbb{N}^*}$, and I remember my teacher confirmed that they did not involve any nasty summation or product.
I remember I was quite out of ideas at the time, and even now I don't really know how to proceed.
Now of course, coupled recursived sequences with constants coefficients such as
$$\left|\begin{array}{lll}\displaystyle a_n=\alpha a_{n-1}+\beta b_{n-1}\\\displaystyle b_n=\gamma a_{n-1}+\delta b_{n-1}\end{array}\right.$$
are swiftly resolved through some linear algebra, provided that the matrix $\begin{pmatrix}\alpha & \beta \\ \gamma & \delta\end{pmatrix}$ is diagonalizable.
But here, the coefficients are not constant. I tried to diagonalize the matrix nonetheless, but even though for any fixed $n\in\mathbb{N}$ it is still diagonalizable, the change of basis matrix (toward the eigenvectors basis) isn't constant with respect to $n$, either. So the straightforward diagonalization method is a dead end.
There might be a purely analytical proof, I don't know.
Any suggestions ?
 A: Hint.
Make 
$$
B_n = \frac{n+2}{(n+1)^2+1}b_n
$$
NOTE
$$
\left(
\begin{array}{c}
a_n\\
B_n
\end{array}
\right) = \left(
\begin{array}{cc}
n & 1\\
-1& -n
\end{array}
\right)\left(
\begin{array}{c}
a_{n-1}\\
B_{n-1}
\end{array}
\right)
$$
Calling now 
$$
M_n = \left(
\begin{array}{cc}
n & 1\\
-1& -n
\end{array}
\right)
$$
for $n$ even we have the curious behavior 
$$
\prod_{k=0}^n M_k = \frac 12 n!\left(
\begin{array}{cc}
-1 & 1\\
-1&  1
\end{array}
\right)
$$
hence for $n$ even
$$
\left(
\begin{array}{c}
a_n\\
B_n
\end{array}
\right) = \frac 12 n!\left(
\begin{array}{cc}
-1 & 1\\
-1&  1
\end{array}
\right)\left(
\begin{array}{c}
a_0\\
B_0
\end{array}
\right)
$$
for $n$ odd is left to the reader.
A: Just computing the very first terms, it is clear that $$a_n=(-1)^n (n-1)!$$
Now, using Cesareo's hint and what you wrote in comments, we then face
$$n B_n=(1-n^2)B_{n-1}-(-1)^n (n-1)!\qquad \text{with} \qquad B_1=3$$ Now, let $B_n=(n-1)! \,C_n$ to get
$$n \,C(n)=-(n+1)\, C(n-1)-(-1)^n\qquad \text{with} \qquad C_1=3$$ to get
$$C_n=(-1)^{n+1}(2n+1)$$
Back to $b_n$,this would give
$$b_n=(-1)^{n+1}\frac{ (2 n+1) (n^2+2n+2) (n-1)!}{n+2}$$
$$\frac{b_n}{n^2a_n}=-\frac{(2 n+1) (n^2+2n+2)}{n^2 (n+2)}=-2-\frac{1}{n}+O\left(\frac{1}{n^2}\right)$$
