might be an over-complicated solution let $a_1+a_2+...+a_n = X_nn^2$ for all $n \geq 2016$ so , we have that $a_n = X_nn^2 - X_{n-1}(n-1)^2 $ for all $n \geq 2016$ from the problem we get the inequality $0 < X_nn^2 - X_{n-1}(n-1)^2 \leq (n+2016)^2 $ for all $n \geq 2016$ let $\lambda = \frac{(2016)^2+2017}{2} $ , first observation is that $2\lambda^2 > (\lambda+2016)^2 $ hence $2n^2 > (n+2016)^2 $ for all $n > \lambda$ so we have $X_n - X_{n-1} \leq 1 $ for all $n > \lambda$ otherwise it will contradict with $RHS$ of the inequality $lemma 1$ , if $ X_n - X_{n-1} = 1$ then $X_{n-1} \leq 2016$ for all $n > \lambda$ : suppose that $X_{n-1} \geq 2017$ then we have $(n+2016)^2 = n^2 + 4032n + 2016^2 \geq X_nn^2 - X_{n-1}(n-1)^2 = n^2 + X_{n-1}(2n-1) \geq n^2 +4034n -2017 $ so $n \leq \frac{(2016)^2+2017}{2} = \lambda$ contradiction ! $lemma 2$, if $X_n - X_{n-1} < 0 $ then $X_{n-1} > \frac{n-2}{2}$ for all $n > \lambda$ since $ X_nn^2 - X_{n-1}(n-1)^2 > 0 $ , $X_nn^2 > X_{n-1}(n-1)^2 \geq (X_n+1)(n-1)^2$ hence $X_n(2n-1) > (n-1)^2 , X_n > \frac{(n-1)^2}{2n-1} > \frac{(n-1)^2}{2n} > \frac{n-2}{2} $ thus $X_{n-1} > X_n > \frac{n-2}{2} , X_{n-1} > \frac{n}{2} $ so , summing all the lemma we have that for all $n > \lambda$ if $X_n - X_{n-1} = 1$ then $X_{n-1} \leq 2016 $ , if $X_n - X_{n-1} < 0 $ then $x_n > \frac{n}{2}$ so we consider three intervals $(-\infty,2016] , (2016,\frac{\lambda}{2}] , (\frac{\lambda}{2},\infty) $ and look at the series $X_{\lambda+1} , X_{\lambda+2} , .....$ if $X_i$ lies in the middle interval , it is easily seen that all $X's$ after ward is equal to $X_i$ if $X_i$ lies in the left interval , then for the next value of $X's$ , it must be non-decreasing and eventually will reach constant since $X's$ are integers if $X_i$ lies in the right interval , then for the next value of $X's$ ,it must be non-increasing which afterward the series would ended up in the first interval or middle interval or ... constant ... which anyway will eventually reach constant hence we could deducted that series $X_i$ will eventually constant thus for large enough $n$ , there exist constant $X$ ,such $a_n = X(2n-1)$ so , $a_n - a_{n-1} = 2X$ for all $n$ large enough which is constant ...

Nice problem. For all $n\in \mathbb{Z}^+$, let $t_n=\frac{a_1+a_2+...+a_n}{n^2}$. Clearly $a_{n+1}=(n+1)^2t_{n+1}-n^2t_n$. The main idea is that $a_{n+1}-a_n$ is eventually constant iff $\{ t_n\}$ is eventually constant. First we will show that there exist finitely many distinct value of $t_k$ such that $t_{k+1} \geq t_k+1$ where $k\in \mathbb{Z}^+$. Suppose to the contary and consider infinite $k>2016$ such that $t_{k+1} \geq t_k+1$. We get $a_{k+1} \geq (k+1)^2(t_k+1)-k^2t_k=(2k+1)t_k+k^2+2k+1$. On the other hand, we have $k+1>2016\implies (k+1+2016)^2 \geq a_{k+1}$. So, $$(k+2017)^2 \geq (2k+1)t_k+k^2+2k+1\implies 4032k+(2017^2-1) \geq (2k+1)t_k.$$And since $(2017^2-1)(2k+1) >4032k+(2017^2-1)$, we get $2017^2-1 >t_k$. This contradict the hypothesis that there exist infinite distinct $k\in \mathbb{Z}^+$ that $t_{k+1} \geq t_k+1$. So, there exist finite set $\mathcal{H}=\{ u\mid t_{u+1} \geq t_u+1 \}$. Then, we will show that there exist constant $C$ such that $t_{\ell+1} \geq t_{\ell}$ for all $\ell >C$ For each positive integer $\ell >2016$ that $t_{\ell+1} <t_{\ell}$, we have $t_{\ell} \geq t_{\ell+1}+1$. Since $a_{\ell}=(\ell+1)^2t_{\ell+1}-\ell^2t_{\ell} >0$, we get $(\ell+1)^2t_{\ell+1} >\ell^2t_{\ell} \geq \ell^2(t_{\ell+1}+1)$. So, \begin{align*} & (2\ell+1)t_{\ell+1} > \ell^2\implies (2\ell+1)(a_1+a_2+...+a_{\ell+1}) >\ell^2(\ell+1)^2\\ & \implies \frac{\ell^2(\ell+1)^2}{2\ell+1} <a_1+a_2+...+a_{\ell+1} \leq (a_1+a_2+...+a_{2015})+\sum_{i=2016}^{\ell+1}{(i+2016)^2}\\ & \implies \ell^2(\ell+1)^2 <(2\ell+1)X+\frac{(2\ell+1)(\ell+1)(\ell+2)(2\ell+3)}{6} +4032\times \frac{(\ell+1)(\ell+2)(2\ell+1)}{2} +2016^2(\ell+1)(2\ell+1)\\ \end{align*}where $X=a_1+a_2+...+a_{2015}-\sum_{i=1}^{2015}{(i+2016)^2}$ is a constant value. Since LHS is a polynomial with degree four with leading coefficient $1$ while RHS is a polynomial degree four with leading coefficient $\frac{2}{3}$, we get that there exists constant $C$ such that when $\ell >C$ the inequality is false. Hence, we finish this part here. From the above two parts, we get that $t_{C+1} \leq t_{C+2} \leq ...$. We, then, have two cases before finishing the problem: There exist infinitely many positive integer $k$ such that $t_{C+k} <t_{C+k+1}$ Let these $k$s are $k_1,k_2,...$. Easy to see that this must be a strictly increasing sequence. But the first part gives us $k_i\in \mathcal{H}$ for all $i$, this is impossible because $\mathcal{H}$ is a finite set. There exist only finitely many positive integer $k$ that $t_{C+k} <t_{C+k+1}$. This means there exists positive integer $D$ such that $t_D=t_{D+1}=...$, so $\{ t_n\}$ is eventually constant, done.

Denote with $k_n=\dfrac{a_1+...+a_n}{n^2}$ for any $n\in \mathbb{N}$. From $(n+1)^2|n^2k_n+a_{n+1}$ we infer $(n+1)^2|(2n+1)k_n-a_{n+1}\ \ (*)$. If $(2n+1)k_n>a_{n+1}$, we infer that $(n+1)^2\le (2n+1)k_n-a_{n+1}<(2n+1)k_n$, so $$a_1+....+a_n=n^2k_n>\dfrac{n^2(n+1)^2}{2n+1}$$However, $$a_1+...+a_n\le 1^2+2^2+...+(n+2017)^2=\dfrac{(n+2017)(n+2018)(2n+2035)}{6}$$which implies $6n^2(n+1)^2<(n+2017)(n+2018)(2n+2035)(2n+1)$. This inequality cannot hold for $n$ big enough (say $n\ge N$) as the leading coefficient in the LHS is bigger than the leading coefficient in the RHS. Thus, for $n\ge N$, $a_{n+1}\ge (2n+1)k_n$. From $(*)$ and the fact that $0\le a_{n+1}-(2n+1)k_n<a_n<2(n+1)^2$, we get that $\left ( a_{n+1}-(2n+1)k_n \right ) \in \{0,(n+1)^2\}$. If $a_{n+1}-(2n+1)k_n=(n+1)^2$, we get that $k_{n+1}=k_n+1$ and if $a_{n+1}-(2n+1)k_n=0$, we get that $k_{n+1}=k_n$. But if $a_{n+1}=(2n+1)k_n+(n+1)^2$, as $a_{n+1}\le (n+2017)^2$, we infer that $k_n<2017$. Therefore, for $n$ big enough we have that $(k_n)_{n\ge M}$ becomes constant, so $b_n=2k_n=\mathrm{constant}$.

Update: I've edit my solution in #4 to be more well-formatted and easier to read.