Bulgarian Mathematical Olympiad 1987 The sequence $x_1,x_2,\dots$ is defined by the equalities $x_1=x_2=1$ and $$x_{n+2}=14x_{n+1}-x_n-4, n\geq 1$$
Prove that each number of the given sequence is a perfect square.

I used the standard way to solve recurrence relations and arrived at
$$x_n=\frac{1}{6}\bigg( (2+\sqrt{3})^{n-\frac{3}{2}} + (2-\sqrt{3})^{n-\frac{3}{2}} \bigg)^2$$
which I think should be quite close. Any idea how to proceed?
 A: First note that $(x_n)_{n\ge 2}$ is strictly increasing, as can be seen from
$$
x_{n+2} = 9 x_{n+1} + (x_{n+1}-x_n) + 4(x_{n+1}-1) \, .
$$
Now – inspired by How do I prove that the recurrence contains no perfect square? – we consider the quadratic form
$$
 Q(x, y) = x^2 - 14xy + y^2 + 4x + 4y + 4 \, .
$$
$Q$ is chosen such that $ Q(x_{n+1}, x_n)$ an invariant: 
$$
 Q(x_{n+2}, x_{n+1}) = Q(14x_{n+1}-x_n-4, x_{n+1}) = Q(x_{n+1}, x_n) \, ,
$$
and since $Q(x_2, x_1) = Q(1,1) = 0$ we have
$$
 Q(x_{n+1}, x_n) = 0 
$$
for $n = 1, 2, 3, \ldots$ . Then also
$$
 Q(x_{n+1}, x_{n+2}) = 0 
$$
because of the symmetry, which means that $y=x_n$ and $y=x_{n+2}$ are (distinct) solutions of the quadratic equation
$$
 0 = Q(x_{n+1}, y) = y^2 + (4-14 x_{n+1})y + (x_{n+1}+2)^2 \, .
$$
Using Vieta's formulas it follows that
$$
 x_n \cdot x_{n+2}  = (x_{n+1}+2)^2
$$
for all $n$, so that $x_{n+2}$ is a perfect if $x_n$ is. And since $x_1$ and $x_2$ are perfect squares, it follows that all $x_n$ are perfect squares.

Remark: $x_n = 1, 1, 9, 121, 1681, 23409, 326041, 4541161, \ldots$ is (apart from the initial element) the sequence A046184 in the  On-Line Encyclopedia of Integer Sequences®:

A046184       Indices of octagonal numbers which are also square. 

A: The equation is equivalent to
$$
x_{n+2}-14x_{n+1}+x_n=-4\tag1
$$
which is the same as
$$
\left(x_{n+2}-\tfrac13\right)-14\left(x_{n+1}-\tfrac13\right)+\left(x_n-\tfrac13\right)=0\tag2
$$
Using the standard methods for Linear Difference Equations, we note that
$$
x^2-14x+1=0\tag3
$$
has roots $7+4\sqrt3$ and $7-4\sqrt3$, so that $(2)$ has solutions
$$
x_n=\tfrac13+c_1\left(7+4\sqrt3\right)^{n-1}+c_2\left(7-4\sqrt3\right)^{n-1}\tag4
$$
Since $x_1=x_2=1$, we get
$$
\begin{align}
x_n
&=\tfrac13+\tfrac{2-\sqrt3}6\left(7+4\sqrt3\right)^{n-1}+\tfrac{2+\sqrt3}6\left(7-4\sqrt3\right)^{n-1}\\[3pt]
&=\tfrac13+\tfrac{2-\sqrt3}6\left(2+\sqrt3\right)^{2n-2}+\tfrac{2+\sqrt3}6\left(2-\sqrt3\right)^{2n-2}\\
&=\left(\tfrac{\sqrt3-1}{2\sqrt3}\left(2+\sqrt3\right)^{n-1}+\tfrac{\sqrt3+1}{2\sqrt3}\left(2-\sqrt3\right)^{n-1}\right)^2\\[3pt]
&=u_n^2\tag5
\end{align}
$$
Since $2+\sqrt3$ and $2-\sqrt3$ are roots of
$$
u^2-4u+1=0\tag6
$$
we have $u_1=u_2=1$ and
$$
u_n=4u_{n-1}-u_{n-2}\tag7
$$
Therefore, $u_n\in\mathbb{Z}$ and $x_n=u_n^2$.

For example, using $(7)$
$$
\{u_n\}=\{1,1,3,11,41,153,571,\dots\}\tag8
$$
and using $(1)$
$$
\{x_n\}=\{1,1,9,121,1681,23409,326041,\dots\}\tag9
$$
A: Long, but hopefully useful version ... Solving recurrence using characteristic polynomials (e.g. like here)
$$\color{red}{x_n}=\frac{1}{6}\left((26-15\sqrt{3})\cdot(7+4\sqrt{3})^n + (26+15\sqrt{3})\cdot(7-4\sqrt{3})^n + 2\right)=\\
\frac{1}{6}\left((26-15\sqrt{3})(2+\sqrt{3})^{2n} + (26+15\sqrt{3})(2-\sqrt{3})^{2n} + 2\right)=\\
\frac{1}{6}\left((2-\sqrt{3})(2+\sqrt{3})^{2(n-1)} + (2+\sqrt{3})(2-\sqrt{3})^{2(n-1)} + 2\right)=\\
\frac{1}{6}\left(\sqrt{2-\sqrt{3}}\cdot(2+\sqrt{3})^{n-1} + \sqrt{2+\sqrt{3}}\cdot(2-\sqrt{3})^{n-1}\right)^2=\\
\frac{1}{6}\left(\frac{\sqrt{3}-1}{\sqrt{2}}\cdot(2+\sqrt{3})^{n-1} + \frac{\sqrt{3}+1}{\sqrt{2}}\cdot(2-\sqrt{3})^{n-1}\right)^2=\\
\frac{1}{6}\left(\frac{3-\sqrt{3}}{\sqrt{6}}\cdot(2+\sqrt{3})^{n-1} + \frac{3+\sqrt{3}}{\sqrt{6}}\cdot(2-\sqrt{3})^{n-1}\right)^2=\\
\left(\frac{3-\sqrt{3}}{6}\cdot(2+\sqrt{3})^{n-1} + \frac{3+\sqrt{3}}{6}\cdot(2-\sqrt{3})^{n-1}\right)^2=\color{red}{(z_{n-1})^2}$$
where
$$z_{n+2}=4z_{n+1}-z_n, \space z_1=1, z_2=3 \tag{1}$$
Note that
$$\left(3-\sqrt{3}\right)\left(2+\sqrt{3}\right)=3+\sqrt{3}$$
and the recurrence in $(1)$ isn't really guessing, since both $2+\sqrt{3}$ and $2-\sqrt{3}$ are roots of $x^2-4x+1=0$, which is the  characteristic polynomial of $(1)$.
