Post here my comment as an answer (actually I'm not sure this is the right answer). Define table of `x` values (just as given by OP): xvals = ParallelTable[ If[TrueQ[ Length[Select[ Select[PowersRepresentations[n, 3, 2], Times @@ # != 0 &], Length[#] == Length[Union[#]] &]] >= 6], n, Nothing], {n, 0, 10000}]; Now we can use this to substitute for `x` when solving the given system of equations, and to speed-up the code we can use `ParallelMap`: ParallelMap[ Solve[{# == n1^2 + n2^2 + n3^2, # == n4^2 + n5^2 + n6^2, # == n7^2 + n8^2 + n9^2, # == n1^2 + n4^2 + n7^2, # == n2^2 + n5^2 + n8^2, # == n3^2 + n6^2 + n9^2}, {n1, n2, n3, n4, n5, n6, n7, n8, n9}, PositiveIntegers] &, xvals[[;; 50]]]; This gives the solution in terms of positive integer numbers, and because of time and memory consumption one can use several ranges of `xvals` (here the first 50 values are used). One can also solve for particular `x`: x = 3051; sol=Solve[{x == n1^2 + n2^2 + n3^2, x == n4^2 + n5^2 + n6^2, x == n7^2 + n8^2 + n9^2, x == n1^2 + n4^2 + n7^2, x == n2^2 + n5^2 + n8^2, x == n3^2 + n6^2 + n9^2}, {n1, n2, n3, n4, n5, n6, n7, n8, n9}, PositiveIntegers] This gives many solutions and `{n1 -> 29, n2 -> 41, n3 -> 23, n4 -> 1, n5 -> 37, n6 -> 41, n7 -> 47, n8 -> 1, n9 -> 29}` is among them. To check if all the solutions are unique we can do: Select[DeleteDuplicates /@ Values[sol], Length[#] == 9 &] In the case of `x = 3051` we have repeated values (1, 29, 41) and this `Select` command gives empty list as an output `{}`, that means there is no unrepeated values.