Does $(n+1)(n-2)x_{n+1}=n(n^2-n-1)x_n-(n-1)^3x_{n-1}$ with $x_2=x_3=1$ define a sequence that is integral at prime indices?
The given difference equation can be solved in the following way. We have for $n\ge 3$, $$ (n-2)(y_{n+1}-y_n) = (n-1)^2(y_n-y_{n-1}),\\ \frac{y_{n+1}-y_n}{n-1}=(n-1)\frac{y_{n}-y_{n-1}}{n-2}. $$ If we let $\displaystyle z_n=\frac{y_{n}-y_{n-1}}{n-2}$, it follows that $$ z_{n+1}=(n-1)z_n=(n-1)(n-2)z_{n-1}=\cdots =(n-1)!z_3=(n-1)!\frac{3x_3-2x_2}{1}=(n-1)! $$ This gives $$ y_{n+1}-y_n=(n-1)(n-1)!=n!-(n-1)!, $$ hence $y_n = nx_n = (n-1)!+c$ for some $c$. Plugging $n=2$ yields $2=1!+c$, thus we have that $$\displaystyle x_n =\frac{(n-1)!+1}{n}.$$
The $n^{\rm{th}}$ term of the sequence is $\dfrac{(n-1)! + 1}{n}$, which is an integer if and only if $n$ is prime (according to Wilson's Theorem).
Too long for a comment:
Numbers don't make any sense
Actually, they do ! :-) Just take a closer look at the sequence's composite-index denominators, and notice the following:
$$ a_{\color{blue}4}=\frac{7}{\color{blue}4},\quad a_{5}=5,\quad a_{\color{blue}6}=\frac{121}{\color{blue}6},\quad a_{7}=103,\quad a_{\color{blue}8}=\frac{5041}{\color{blue}8},\quad a_{\color{blue}9}=\frac{40321}{\color{blue}9},\quad \\~\\~\\ a_{\color{blue}{10}}=\frac{362881}{\color{blue}{10}}, a_{11}=329891,\quad a_{\color{blue}{12}}=\frac{39916801}{\color{blue}{12}},\quad a_{13}=36846277,\quad\dots $$
Let us now rewrite the sequence's prime-indexed elements in a similar manner, for a more uniform approach:
$$ a_{\color{blue}4}=\frac{7}{\color{blue}4},\quad a_{\color{blue}5}=\frac{25}{\color{blue}5},\quad a_{\color{blue}6}=\frac{121}{\color{blue}6},\quad a_{\color{blue}7}=\frac{721}{\color{blue}7},\quad a_{\color{blue}8}=\frac{5041}{\color{blue}8},\quad a_{\color{blue}9}=\frac{40321}{\color{blue}9},\quad \\~\\~\\ a_{\color{blue}{10}}=\frac{362881}{\color{blue}{10}}, a_{\color{blue}{11}}=\frac{3628801}{\color{blue}{11}},\quad a_{\color{blue}{12}}=\frac{39916801}{\color{blue}{12}},\quad a_{\color{blue}{13}}=\frac{479001601}{\color{blue}{13}},\quad\dots $$
At this point, we might be able to cheat, and use OEIS to identify the afferent sequence $\color{blue}{b_n=n\cdot a_n}$ by its first few elements, yielding three possible suspects: but let's say that our mathematical virtue and intellectual integrity will prevail over our base urges of rampant curiosity, and we might resist the temptation to do so. Could we then, without any outside aide, still manage to deduce an expression for $b_n$ ? Indeed, even a superficial glance will undoubtedly reveal the growth to resemble what one might otherwise expect to see in a geometric progression, rather than an arithmetic one. We then have:
$$ r_{\color{blue}4}=\frac{25}{7}\simeq\color{blue}4,\quad r_{\color{blue}5}=\frac{121}{25}\simeq\color{blue}5,\quad r_{\color{blue}6}=\frac{721}{121}\simeq\color{blue}6,\quad r_{\color{blue}7}=\frac{5041}{721}\simeq\color{blue}7,\quad\dots $$
In short, $\color{blue}{r_n=\dfrac{b_{n+1}}{b_n}\simeq n}.~$ A type of “modified geometric progression”, as it were, with a variable ratio equal to the index itself: Does this sound in any way familiar ? If no, there is still no need to worry: We'll get to the bottom of it soon enough, anyway. More precisely, we have $\color{blue}{b_{n+1}=n\cdot b_n-(n-1)}.~$ Now we are left with showing that, for prime values of the index p, $~\color{blue}{a_p=\dfrac{b_p}p\in\mathbb N},~$ starting from $\color{blue}{b_2=2}.~$ Could mathematical induction work in this case, or does the seemingly random distribution of primes throw a wrench into such an approach ? And, if so, could we then maybe slightly modify it, to fit the new conditions ? What would you say ? :-)