Four squares such that the difference of any two is a square?

I could manage to find many almost solutions (with the aim to derive full solutions) by a systematic search:

My first trials, code and ideas are posted here, where I received very helpful input that directed me to a promising approach, proposed by Peter: Search for 6-tuples $(s,t,u,t+u,t+u−s,t−s)$ consisting of perfect squares. I implemented a Python script which systematically searched for those 6-tuples and it found a lot of instances. Unfortunatelly at a certain point it became slow and did not scale anymore.

With the help of Arty, who developed an enhaced version, I could manage to find a first allmost matches. This performance breakthrough origins from this question, I asked on SO.

On a machine with 128GB RAM, an Intel Xeon CPU E5-2630 v4, 2.20GHz processor (and two graphic cards of type Tesla V100 with 16GB RAM) I generated a large textfile of these 6-tuples $(s,t,u,t+u,t+u−s,t−s)$.

To make the final step, namely to identify these four squares, I implemented the Mathematica Script pythagorean.nb, which outputs an "almost solution": [w=40579, x=-58565, y=-65221].

From the data set I selected the corresponding row:

42228, 51060, 185472, 1783203984, 2607123600, 34399862784, 37006986384, 35223782400, 823919616

Hence $s=42228^2,t=51060^2,u=185472^2$. To calculate the last variable $z$, I need to use the equation $u+y^2=z^2$ leading to $z=196605,294$.

Now my idea is to multiply all integers $w,x,y,z$ with a suitable integer, for example with $10^{10}$ to bring us closer towards a solution:

• $z^2-y^2=1854720000103154^2$
• $z^2-x^2=1876800000101940.2^2$
• $z^2-w^2=1923720000099454^2$
• $y^2-x^2=287040000000000^2$
• $y^2-w^2=510600000000000^2$
• $x^2-w^2=422280000000000^2$

An excerpt of the list containing such "almost" solutions, I could generate so far is (raw output of Mathematica):

• $w=40579,x=-58565,y=-65221,z=196605.2940$
• $w=53900,x=-71540,y=-91756,z=157415.4376$
• $w=534240,y=-540600,z=-554115,x=537302.7941$
• $w=67375,x=-89425,y=-114695,z=196769.2970$
• $w=81158,x=-117130,y=-130442,z=393210.5880$

The Mathematica Code is shown below:

arr = Import[
   "C:/Users/esultano/git/pythagorean/pythagorean_stu_Arty_.txt", 
   "CSV", "HeaderLines" -> 0];
f[i_] := Part[arr[[i]], 1 ;; 3];
len = Length[arr];
For[i = 1, i < len, i++, {
  triplet = f[i];
  s = triplet[[1]];
  t = triplet[[2]];
  u = triplet[[3]];
  ins = FindInstance[
    x*x - w*w == s*s && y*y - w*w == t*t && w != 0, {w, x, y}, 
    Integers];
  If[Length[ins] > 0, 
   Print[Append[First[ins], 
     z -> N[Sqrt[u^2 + First[ins][[3, 2]]^2], 10]]], Continue];
  ins = FindInstance[
    y*y - w*w == t*t && z*z - y*y == u*u && w != 0, {w, y, z}, 
    Integers];
  If[Length[ins] > 0, 
   Print[Append[First[ins], 
     x -> N[Sqrt[s^2 + First[ins][[1, 2]]^2], 10]]], Continue];
  }
 ]