When the $a^2+b+c+d,b^2+a+c+d,c^2+a+b+d,d^2+a+b+c$ are all perfect squares? Find all $a,b,c,d\in \mathbb{Z}^+$, which $a^2+b+c+d,b^2+a+c+d,c^2+a+b+d,d^2+a+b+c$ are all perfect squares.
I found $(1,1,1,1)$, but I can't find more.
Is $a=b=c=d$ true?
 A: Since the expressions are symmetric, then we can set $a\leq b\leq c\leq d$ without loss of generality. There is an infinite number of positive integer solutions if we assume $a=1$ and $b=c=d=n$ which yields,
$$\begin{aligned}
a^2+b+c+d &= 3n+1\\
a+b^2+c+d &= (n+1)^2\\
a+b+c^2+d &= (n+1)^2\\
a+b+c+d^2 &= (n+1)^2
\end{aligned}\tag1$$
and it is easy to find solve $3n+1 = y^2$. 
Excluding this infinite family, it seems the only solutions with all variables $1\leq a\leq b\leq c\leq d < 100$ are,
$$a,b,c,d = 6,6,11,11$$
$$a,b,c,d = 40,57,96,96$$
though I am not sure of this result. (Anyone can verify it?)
A: The solutions are $(11, 11, 6, 6)$, $(96, 96, 57, 40)$, and $(x, x, x, 1)$, where $3x + 1$ is a square.
Without loss of generality, $a \ge b \ge c \ge d$, since it is symmetric (not just cyclic). Then$$a^2 < a^2 + b + c + d \le a^2 + 4a < (a+2)^2 \implies a^2 + b + c + d = (a+1)^2.$$Thus, $b + c + d = 2a + 1$. Then$$b \ge {1\over3}(2a + 1) \implies a \le {3\over2}b,$$so$$b^2 < b^2 + c + d + a \le b^2 + {7\over2}b < (b+2)^2 \implies b^2 + c + d + a = (b+1)^2.$$Thus, $b + c + d = 2a + 1$ and $c + d + a = 2b + 1$, which implies $a = b$, thus $c + d = a + 1$.
Consequently,$$c \ge{1\over2}(2a+1) \implies a < 2c,$$so$$(c+1)^2 \le c^2 + 2a + d \le c^2 + 5c < (c+3)^2,$$and hence we have two cases.
Case 1. If $c^2 + 2a + d = (c+1)^2$, we must have $a = b = c$ and $d = 1$. Thus, we find that $(a, b, c, d) = (x, x, x, 1)$ where $x$ is such that $3x + 1$ is a square. This is one set of solutions.
Case 2. If $c^2 + 2a + d = (c + 2)^2$, so $2a + d = 3c + 4$. Combined with $c + d = 2a + 1$, we derive that $a = (5/2)d - 4$ and $c = (3/2)d - 3$. So$$(d+1)^2 \le d^2 + 2a + c = d^2 + {{13}\over2}d - 11 < (d+4)^2.$$So setting this equal to $(d+1)^2$, $(d+2)^2$, $(d+3)^2$, we find the integer solutions $d = 6$ and $d = 40$. This gives the solutions $(11, 11, 6, 6)$ and $(96, 96, 57, 40)$.
Hence, the solutions are $(11, 11, 6, 6)$, $(96, 96, 57, 40)$, and $(x, x, x, 1)$, where $3x + 1$ is a square.
A: I ran a small program in C# to use brute force to find many answers:
1,5,5,5
1,8,8,8
5,1,5,5
5,5,1,5
5,5,5,1
8,1,8,8
8,8,1,8
8,8,8,1
96,96,57,40
Here's the program:
using System;

namespace PerfectSquare
{
    class MainClass {
        public static void Main (string[] args) {
            long a,b,c,d;
            const int range = 100;  
            for (a=1; a<=range; a++) {
                for(b=1;b<=range;b++) {
                    for(c=1;c<=range;c++) {
                        for(d=1;d<=range;d++) {
                            if (!IsPerfectSquare(a*a+b+c+d)) continue;
                            if (!IsPerfectSquare(b*b+a+c+d)) continue;
                            if (!IsPerfectSquare(c*c+a+b+d)) continue;
                            if (!IsPerfectSquare(d*d+a+b+c)) continue;
                            Console.WriteLine (String.Format ("{0} {1} {2} {3}",a,b,c,d));
                        }
                    }
                }
            }
            Console.WriteLine ("Finished.");
        }
        private static bool IsPerfectSquare(long value) {
            System.Double test;
            test = Math.Sqrt(value);
            return Math.Floor(test) == test;
        }
    }
}

A: Do brute force search. Here is a piece of Python code that should do the job
from math import sqrt, floor

_is_whole = lambda x: floor(x) == x
_is_ps = lambda n: _is_whole(sqrt(n))

min_int = 1
max_int = 150
for a in range(min_int, max_int + 1):
    for b in range(a, max_int + 1):
        for c in range(b, max_int + 1):
            for d in range(c, max_int + 1):
                for x in [a ** 2 + b + c + d, a + b ** 2 + c + d,
                          a + b + c ** 2 + d, a + b + c + d ** 2]:
                    if not _is_ps(x):
                        break
                else:
                    print a, b, c, d

When run, it should produce the following output:
1 1 1 1
1 5 5 5
1 8 8 8
1 16 16 16
1 21 21 21
1 33 33 33
1 40 40 40
1 56 56 56
1 65 65 65
1 85 85 85
1 96 96 96
1 120 120 120
1 133 133 133
6 6 11 11
40 57 96 96
