Let $b_n = n^2a_n + 2$ for all $n \ge 1$. The recurrence relation becomes
$$(n + 1)(b_{n+1} - 2) = 2(2n + 1)(b_n - 2) + 2(3n + 1).$$
This is the same as
$$b_{n+1} = \frac{2(2n + 1)}{n + 1}b_n.$$
It follows that
$$b_n = \frac{2(2n-1)}{n} b_{n-1} = \frac{2^2(2n-1)(2n-3)}{n(n-1)} b_{n-2} = \dots = \frac{2^{n-1}(2n-1)!!}{n!} b_1,$$
where $(2n-1)!! = (2n-1)(2n-3)\cdots(1)$. Since $b_1 = a_1 + 2 = 2$, we have
$$b_n = \frac{2^n(2n-1)!!}{n!} = \frac{(2n)!!(2n-1)!!}{(n!)^2} = \frac{(2n)!}{(n!)^2} = \binom{2n}{n}.$$
Therefore, we obtain
$$a_n = \frac{b_n - 2}{n^2} = \frac{1}{n^2} \left[ \binom{2n}{n} - 2 \right].$$
By Vandermonde's identity, we can rewrite this as
$$a_n = \frac{1}{n^2} \left( \sum_{k=0}^{n} \binom{n}{k}^2 - 2 \right) = \frac{1}{n^2} \left( \sum_{k=1}^{n-1} \binom{n}{k}^2 \right).$$

When $n = p$ is a prime, we have $p \mid \binom{p}{k}$ for $1 \le k \le n-1$, and so $p^2 \mid \binom{p}{k}^2$. This implies $a_p \in \mathbb{Z}$. Clearly, $a_p > 0$. Since there are infinitely many prime numbers, there are infinitely many terms which are positive integers.