The quick version is $n_0 = 0, \; \; n_1 = 40,$ then
$$ \color{magenta}{ n_{k+2} = 98 n_{k+1} - n_k + 40}. $$
Given an $(x,y)$ pair with $3x^2 - 2 y^2 = 1$ we then take $n = (x^2-1)/ 2 = (y^2 - 1)/ 3. $
The first few $x,y$ pairs are
$$ x=1, \; y= 1 , \; n=0 $$
$$ x=9, \; y=11, \; n=40 $$
$$ x= 89, \; y=109, \; n=3960 $$
$$ x=881, \; y=1079, \; n= 388080 $$
$$ x=8721, \; y=10681, \; n= 38027920 $$
$$ x=86329, \; y=105731, \; n= 3726348120 $$
and these continue forever with
$$ x_{k+2} = 10 x_{k+1} - x_k, $$
$$ y_{k+2} = 10 y_{k+1} - y_k. $$
$$ n_{k+2} = 98 n_{k+1} - n_k + 40. $$
People seem to like these recurrences in one variable. The underlying two-variable recurrence in the pair $(x,y)$ can be abbreviated as
$$ (x,y) \; \; \rightarrow \; \; (5x+4y,6x+5y) $$ beginning with
$$ (x,y) = (1,1) $$
The two-term recurrences for $x$ and $y$ are just Cayley-Hamilton applied to the matrix
$$ A \; = \;
\left( \begin{array}{rr}
5 & 4 \\
6 & 5
\end{array}
\right) ,
$$
that being
$$ A^2 - 10 A + I = 0. $$
You may as well assume that the first term is a square. Suppose that
$$M, M + aD, M + bD, M + cD, M + dD$$
are equal to $v^2$, $w^2$, $x^2$, $y^2$, and $z^2$. Then one obtains a solution to the equations
$$\frac{w^2 - v^2}{a} = \frac{x^2 - v^2}{b} = \frac{y^2 - v^2}{c} = \frac{z^2 - v^2}{d}.$$
We are asking whether this has any solutions in integers $w,v,x,y,z$. Thinking of $[w;v;x;y;z] \in \mathbf{P}^4$, the equations above cut out a curve $C$ which is the intersection of three quadrics (the pairwise differences of the first term and the second, third, and fourth term). Assuming $a,b,c,d$ are distinct, the curve $C$ is smooth of has genus $5$, and hence has only finitely many rational points by Faltings' theorem. Hence there will only be finitely many arithmetic progressions (up to scaling) with at least $5$ squares.
On the other hand, I claim there are infinitely many arithmetic progressions (up to scaling) for which $M$, $M+D$, $M+2D$, and $M+4D$ are all squares. In this case, we obtain the equations:
$$w^2 - v^2 = \frac{x^2 - v^2}{2} = \frac{y^2 - v^2}{4}.$$
This is the intersection $X$ of two quadrics in $\mathbf{P}^3$, which defines a curve of genus one. The curve corresponding to $M$, $M+D$, $M+2D$, and $M+3D$ all being squares turns out not to have rational points (Fermat) but this one is an elliptic curve of rank one.
In fact, $X$ is (clearly) isogenous to
$$Y^2 = (1+X)(1+2X)(1+4X),$$
which is an elliptic curve of conductor $192$. If $P = [0,1]$, then the odd multiples of $P$ give rise to such sequences, e.g.:
$$P: M = 1, \quad D = 0,$$
$$3P: M = 49, \quad D = 120,$$
$$5P: M = 17497489, \quad D = -4269720,$$
$$7P: M = 4085374755361, \quad D = 82153503191760,$$
etc.
Summary: There will be infinitely many arithmetic progressions (even of length $5$) with $4$ squares, however, for any fixed length (say $10$, or $100$, or $1000$) there will only be finitely many arithmetic progressions up to scaling with at least $5$ squares.
To find the finitely many exceptions with $5$ or more squares could be a little bit annoying. For example, one could find all the rational points on the $\binom{9}{4}$ genus $5$ curves coming from assuming the first term and four more are all squares. Some of these will be duplicates by symmetry. It might first make sense to consider the $\binom{9}{3}$ genus one curves coming from assuming the first term and three more are all squares, see which ones have positive rank (although $P = [0,1]$ is a point on $Y^2 = (1+aX)(1+bX)(1+cX)$ and may well have infinite order unless $[a;b;c] = [1;2;3]$). If Michael Stoll is still around, he could probably come up with a complete answer. As far as I know, there may not be any with $5$ squares, although clearly this is false if one replaces $10$ (say) by $25$, since then $1,2,3,4,\ldots,25$ has $5$ squares. (This example also had me confused as to why you had a hard time finding an arithmetic progression with three squares, since $1,2,3,4,5,6,7,8,9,10$ does the job.)
Best Answer
First question is a well known problem studied by Euler and some variants of it still out of our reach . For every positive integers $(x,y,z)$ let $P(x,y,z)$ denote the predicate: $(xy+1)(yz+1)(zx+1)$ is a square and not all of $(xy+1)$, $(yz+1)$,$(zx+1)$ are squares
Proof by infinite descent:
we well prove by the infinite descent that if there exists a solution $(x,y,z)$ such that $P(x,y,z)$ is true than there exists another solution $(p,q,r)$ such that $P(p,q,r)$ is true and $p+q+r$ less than $x+y+z$.
Let $(x,y,z)$ be tuple of positive integers such that $P(x,y,z)$ is true and $x\leq y \leq z$, consider the two positive integers $s_{+}$ and $s_{-}$ defined by : $$s_{\mp}=x+y+z+2xyz\mp\sqrt{(xy+1)(yz+1)(zx+1)} \tag 1$$ This integers verify : $$x^2+y^2+z^2+s_{\mp}^2-2(xy+yz+zx+xs_{\mp}+ys_{\mp}+zs_{\mp})-4xyzs_{\mp}-4=0 \tag2$$ And we can also prove the following important identities (they are basically the same): $$ \begin{align} (x+y-z-s_{\mp})^2&&=&& 4(xy+1)(zs_{\mp}+1) \tag 3 \\ (x+z-y-s_{\mp})^2&&=&& 4(xz+1)(ys_{\mp}+1) \tag 4 \\ (x+s_{\mp}-z-y)^2&&=&& 4(yz+1)(xs_{\mp}+1) \tag 5 \end{align}$$ we can use this identities to proof that $P(x,y,s_{\mp})$ holds, by multiplying $(4)$ and $(5)$ you will get that $(xy+1)(ys_{\mp}+1)(xs_{\mp}+1) $ is a square and not all of $(xy+1)$,$(xs_{\mp}+1)$ , $(ys_{\mp}+1)$ are square (as we have $(xz+1)$ is a square iff $xs_{\mp}+1)$ is a square and $(yz+1)$ is a square iff $(ys_{\mp}+1)$ is a square using $(4)$ and $(5)$).
Now the important par is to prove that either $x+y+s_{+}<x+y+z$ or $x+y+s_{-}<x+y+z $ this is equivalent tp proving that $s_{-}s_{+}<z^2$ which is true because : $$s_{-}s_{+}=x^2+y^2+z^2-2(xy+yz)-4=z^2-x(2z-x)-y(2z-y)-4<z^2 $$ (remember that $x\leq y \leq z$)
Reference: When Is (xy + 1)(yz + 1)(zx + 1) a Square? Kiran S. Kedlaya Mathematics Magazine Vol. 71, No. 1 (Feb., 1998), pp. 61-63
Second Question: The following Pell equation have infinitely many solutions : $$x^2-3y^2=1 $$
For any solution $(n,m)$ of this Pell equation, let $(a,b,c)=(2n-m,2n,2n+m)$ then $ab+1,bc+1,ca+1$ are all squares.