Here are some recipes for finding infinitely many solutions.
The initial setup
It is equivalent to find two nontrivially different arithmetic progressions of three squares, $u^2, v^2, w^2$ and $U^2, V^2, W^2$, where $u,v,w,U,V,W$ are rational, with the same common difference. Once we do, we can substitute $$p = \frac{U+u}{2},\; q = \frac{U-u}{2},\; a = \frac{V+v}{2},\; b = \frac{V-v}{2},\; r = \frac{W+w}{2},\; s = \frac{W-w}{2}$$ and obtain a rational solution to the Diophantine equation we started with, because $(p+q)^2, (a+b)^2, (r+s)^2$ will be $U^2, V^2, W^2$ while $(p-q)^2, (a-b)^2, (r-s)^2$ will be $u^2, v^2, w^2$. We can then scale it up to clear denominators and get integer solutions.
(We can also always change signs of $p,q,a,b,r,s$ to make them all positive, if necessary, so don't worry about that part. First, observe that we can multiply both numbers in a pair by $-1$ and change nothing, so we can assume $p,a,r$ are all positive. Because $pq=ab=rs$, we know $q,b,s$ are all positive or all negative; in the latter case, negate $q,b,s$ and we simply swap our two arithmetic progressions.)
Given that this is what we're thinking about, we naturally begin to re-read Keith Conrad's article on Arithmetic Progressions of Three Squares.
Using elliptic curves
Note: this method is more complicated; however, there is no "waiting time" involved. You go through all the calculations exactly once, and get a solution; you can keep going on the elliptic curve, and generate infinitely many solutions. This has theoretical importance - it also proves that there are in fact infinitely many solutions! However, there's no guarantee that you'll generate all possible solutions in this way.
We see (in Corollary 3.3) that the triples of rational numbers $(u,v,w)$ whose squares are in an arithmetic progression with difference $n$ are in bijection with rational points on the elliptic curve $y^2 = x^3 - n^2 x$ with $y \ne 0$. Also (from Theorem 3.6), all such points have infinite order on the elliptic curve, which means that as soon as we find one triple $(u,v,w)$, we can find infinitely many candidates $(U,V,W)$. The recipe is to
- Convert $(u,v,w)$ it to a point $(x,y)$ on the elliptic curve, using Corollary 3.3.
- Using the group law on elliptic curves, repeatedly add the point to itself any number of times you like.
- Convert the resulting point back to a triple $(U,V,W)$.
This is illustrated in Example 3.8 of the article, which really covers almost everything we need to do in order to find solutions to our problem. The only remaining step is the reduction at the beginning of this answer, which is much easier.
For example, suppose we started with the triple $(u,v,w) = (1,5,7)$, converted it to the point $P = (-12,72)$ on $y^2 = x^3 - 24^2 x$, added it to itself to get $2P = (25, -35)$, and then converted it back to get $(U,V,W) = (\frac{1151}{70}, -\frac{1201}{70}, \frac{1249}{70})$. (All this is directly taken from Example 3.8.) To save time later, we can immediately scale up both triples by $140$, getting $(u,v,w) = (140, 700, 980)$ and $(U,V,W) = (2302, 2402, 2498)$. We then obtain $$(p,q,a,b,r,s) = (1221, 1081, 1551, 851, 1739, 759)$$ by the correspondence in my first paragraph. Indeed, we have an arithmetic progression $$(1221 + 1081)^2, (1551^2 + 851)^2, (1739^2 + 759)^2$$ (with common difference $470400$) and that $1221 \cdot 1081 = 1551 \cdot 851 = 1739 \cdot 759 = 1319901$.
To find the starting $(u,v,w)$ to use here, see the next section.
Using rational points on a circle
Note: this method is less complicated; however, it is also less certain. In practice, we find several solutions relatively quickly, but there is no guarantee of anything. We will eventually find all the solutions, but in the same way as we would eventually find all the solutions by doing a brute-force search. This method is faster than a brute-force search, at least.
Section 2 of Keith Conrad's article also outlines how to find triples of squares in arithmetic progression to begin with. The elliptic curve approach is necessary if you want to be "proactive" about finding solutions $(p,q,a,b,r,s)$. However, if we merely want to get to all of them sooner or later, we can simply generate arithmetic progressions $u^2, v^2, w^2$ and tabulate their common differences until we encounter a repeat. (We might want to replace the common difference by its largest square-free divisor, since scaling $(u,v,w)$ can multiply that common difference by any perfect square.)
In more detail: for every rational number $m \in (1, 1+\sqrt2)$, $\left(\frac{1+2m-m^2}{2}, \frac{m^2+1}{2}, \frac{m^2+2m-1}{2}\right)$ is a positive triple whose squares are in arithmetic progression with common difference $m^3-m$. (The range on $m$ is to avoid duplicates, which are caused by obtaining essentially the same triple where some terms have a different sign.) This generates a multiple of every possible such triple.
If we find two triples generated by rational numbers $m, m'$ such that $\sqrt{\frac{m^3-m}{m'^3-m'}}$ is rational, that means we can scale the triples generated in this way so that the common difference of the APs is the same. This gives us another pair of $(u,v,w)$ and $(U,V,W)$ as in the previous solution.
Going through the rationals between $1$ and $1+\sqrt2$ in order of denominator, the first such pair I find is $m=\frac73$, $m'=\frac75$, which correspond to the triples $(\frac19, \frac{29}{9}, \frac{41}{9})$ and $(\frac{23}{25}, \frac{37}{25}, \frac{47}{25})$; the squares have common difference $\frac{280}{27}$ and $\frac{168}{125}$, respectively. Multiplying the first triple by $9$ and the second by $25$, we get $(u,v,w) = (1,29,41)$ and $(U,V,W)=(23,37,47)$, where $v^2-u^2$, $w^2-v^2$, $V^2-U^2$, and $W^2-V^2$ are all $840$. This, in turn, gives us the solution $(p,q,a,b,r,s) = (12, 11, 33, 4, 44, 3)$, which - as others have already observed - is the first nontrivial solution.
A table
Here are the first few solutions obtained by the method of rational points on a circle, together with the rational numbers $m$ and $m'$ that generate them:
\begin{array}{cc|cccccc}
m & m' & p & q & a & b & r & s \\
\hline
\frac{7}{3} & \frac{7}{5} & 12 & 11 & 33 & 4 & 44 & 3 \\
\frac{7}{3} & \frac{8}{7} & 99 & 95 & 171 & 55 & 209 & 45 \\
\frac{6}{5} & \frac{11}{5} & 105 & 91 & 195 & 49 & 245 & 39 \\
\frac{7}{5} & \frac{8}{7} & 143 & 51 & 187 & 39 & 221 & 33 \\
\frac{7}{6} & \frac{14}{13} & 275 & 62 & 310 & 55 & 341 & 50 \\
\frac{10}{7} & \frac{17}{10} & 209 & 58 & 418 & 29 & 551 & 22 \\
\frac{12}{7} & \frac{19}{16} & 649 & 357 & 1003 & 231 & 1239 & 187 \\
\frac{13}{7} & \frac{13}{8} & 165 & 41 & 451 & 15 & 615 & 11 \\
\frac{13}{7} & \frac{15}{13} & 99 & 68 & 153 & 44 & 187 & 36 \\
\frac{9}{8} & \frac{17}{8} & 341 & 294 & 539 & 186 & 651 & 154 \\
\frac{11}{8} & \frac{19}{8} & 182 & 175 & 490 & 65 & 650 & 49 \\
\frac{13}{8} & \frac{15}{13} & 437 & 231 & 627 & 161 & 759 & 133 \\
\end{array}
A nice subfamily of solutions
If we want a second, more elementary proof that there are infinitely many primitive solutions, we can spot a subsequence inside these examples: the subsequence where $m' = m+1$. This gives us a pair where $\frac{m'^3-m'}{m^3-m}$ is a perfect square if we take $m = \frac{k^2+2}{k^2-1}$ and $m' = \frac{2k^2+1}{k^2-1}$. All integers $k\ge 3$ give us a solution in this way, although not all solutions are of this form.
Going through the whole process and simplifying a little, we get the parameterized family
\begin{align}
p &= (k^2+k+1)(2k+1)(k-2)(k+1) \\
q &= (k^2-k+1)(2k-1)(k+2)(k-1) \\
a &= (k^2-k+1)(2k+1)(k+2)(k+1) \\
b &= (k^2+k+1)(2k-1)(k-2)(k-1) \\
r &= (k^2+k+1)(2k-1)(k+2)(k+1) \\
s &= (k^2-k+1)(2k+1)(k-2)(k-1)
\end{align}
The factored form makes it easy to see that $pq = ab = rs$; some algebra lets us check that $(p+q)^2, (a+b)^2, (r+s)^2$ are in an arithmetic progression with common difference $48 k^2 (k^2 - 1) (k^2 + 2) (2 k^2 + 1)$. The triples generated by this family aren't actually quite primitive; it seems that we can divide through by one of $2$, $4$, $18$, or $36$, depending on $k \bmod 12$. However, no two of them are multiples of each other, so when we do divide through by the GCD, we get infinitely many primitive solutions.
