Limit of recursive sequence $n^2q_n=1+(n-1)^2q_{n-1}+2(n-2)q_{n-2}$

When looking at this riddle, I came across the following sequence for the frequency of sampled integers between 1 and $n$ in a without replacement/without neighbour sampling: $$q_1=1,\quad q_2=1/2,\quad q_n=\frac{1}{n^2}+\frac{(n-1)^2}{n^2}q_{n-1}+\frac{2(n-2)}{n^2}q_{n-2}\quad (n>2)$$ and I wonder if there is a generic mathematical approach to computing an analytic solution for $$\lim_{n\to\infty} q_n$$Assessing the numerical limit with an R code like

q=rep(1,1e7)
for (n in 3:1e7) q[n]=(1+2*q[n-2]+(n-1)*q[n-1])/n
q[1e7]/1e7

led me to $0.432332...$ And a probabilistic reasoning indicates that $(n>1)$ $$\frac{1}{3}\le q_n\le \frac{1}{2}$$

Let $a_n = (n+1)q_{n+1}$. We have $a_0 = a_1 = 1$ and

$$(n+1)a_n = 1 + na_{n-1} + 2a_{n-2}\quad\text{ for }\quad n > 1\tag{*1}$$

Let $f(z) = \sum\limits_{n=0}^\infty a_n z^n$, multiply $(*1)$ by $z^n$ and start to sum from $n = 2$, we get

$$\begin{align} & \left(z\frac{d}{dz} + 1 \right)(f(z) - 1 - z) = \frac{z^2}{1-z} + \left(z\frac{d}{dz}\right)(z(f(z)-1)) + 2z^2f(z)\\ \iff & zf' + f - 1 -2z = \frac{z^2}{1-z} + z^2 f' + zf - z + 2z^2 f\\ \iff & z(1-z)f' + (1-z-2z^2)f = \frac{1}{1-z} \end{align} $$ Solving the ODE give us

$$f(z) = \frac{1-e^{-2z}}{2z(1-z)^2} = \frac{A}{(1-z)^2} + \frac{B}{(1-z)} + g(z)\tag{*2}$$ where $\displaystyle\;\begin{cases} A &= \frac{1-e^{-2}}{2}\\ B &= \frac{1-3e^{-2}}{2}\\ \end{cases}\; $ and $g(z) = \sum\limits_{n=0}^\infty g_n z^n$ is some function analytic over all $\mathbb{C}$.

Expanding $(*2)$ as a power series and compare coefficients of $z^n$ on both sides, we get

$$(n+1)q_{n+1} = a_n = (n+1)A + B + g_n$$

Since $g(z)$ is entire, its coefficients of taylor expansion $g_n$ is bounded. As a result, $$ \begin{align} \lim_{n\to\infty} q_{n+1} &= A + \lim_{n\to\infty} \frac{B + g_n}{n+1} = A = \frac{1 - e^{-2}}{2}\\ &\approx 0.432332358381693654053000252513757798296184227045212059265 \end{align} $$