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? My son gave me the following recurrence formula for $x_n$ ($n\ge2$):

$$(n+1)(n-2)x_{n+1}=n(n^2-n-1)x_n-(n-1)^3x_{n-1}\tag{1}$$
$$x_2=x_3=1$$

The task I got from him:   


*

*The sequence has an interesting property, find it out.

*Make a conjecture and prove it.


Obviously I had to start with a few values and calculating them by hand turned out to be difficult. So I used Mathematica and defined the sequence as follows:
b[n_] := b[n] = n (n^2 - n - 1) a[n] - (n - 1)^3 a[n - 1];
a[n_] := a[n] = b[n - 1]/(n (n - 3));
a[2] = 1;
a[3] = 1;

And I got the following results:
$$ a_4=\frac{7}{4}, \ a_5=5, \ a_6=\frac{121}{6}, \ a_7=103, \ a_8=\frac{5041}{8}, \ a_9=\frac{40321}{9}, \\ a_{10}=\frac{362881}{10}, \ a_{11}=329891, \ a_{12}=\frac{39916801}{12}, \ a_{13}=36846277, \ a_{14}=\frac{6227020801}{14}\dots$$
Numbers don't make any sense but it's strange that the sequence produces integer values from time to time. It's not something that I expected from a pretty complex definition like (1).
So I decided to find the values of $n$ producing integer values of $a_n$. I did an experiment for $2\le n \le 100$:
table = Table[{i, a[i]}, {i, 2, 100}];
integers = Select[table, (IntegerQ[#[[2]]]) &];
itegerIndexes = Map[#[[1]] &, integers]

...and the output was:
{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 
59, 61, 67, 71, 73, 79, 83, 89, 97}

Conjecture (pretty amazing, at least to me):

$a_n$ is an integer if and only if $n$ is prime.

Interesting primality test, isn't it? The trick is to prove that it's correct. I have started with the substitution: 
$$y_n=n x_n$$
...which simplifies (1) a bit:
$$(n-2)y_{n+1}=(n^2-n-1)y_n-(n-1)^2y_{n-1}$$
...but I did not get much further (the next step, I guess, should be rearrangement). 
 A: 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}.$$
A: 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).
A: My process to find the solution probably involves too far a leap of faith, but I will post it anyway.
First, after writing out at least 8 terms, I noticed that, apart from integers, all terms are of the form $a_n=?/n$. So naturally I turned towards the sequence $b_n\equiv a_n\cdot n$, which goes like $3, 7, 25, 121, 721, 5041, 40321, \dots$.
Now my intuition senses the number sequence $504*, 4032*$, which resembles the factorial of 7&8, moreover, a lot of them seems to end in $1$. So I made a further conjecture that  

$$a_n=\frac{(n-1)!+1}{n}.$$

A quick check of Wilson's theorem shows that it is compatible with the original conjecture. Thus, we seem to be on the right track. This can be easily proved by induction.
A: 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}},\quad
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 ? :-)
