(A while ago I reported here doing this in Macsyma.) (Barely noteworthy, but undeserving of oblivion.) Find the 3x3 msqs just using Solve. The square: In[2]:= mat = {{a, b, c}, {d, e, f}, {g, h, i}}; The row and diagonal sums: In[6]:= lin = 15 == # & /@ Flatten[{Total /@ mat, Total /@ Transpose[mat], Tr[mat], Tr[Reverse[mat]]}] Out[6]= {15 == a + b + c, 15 == d + e + f, 15 == g + h + i, 15 == a + d + g, 15 == b + e + h, 15 == c + f + i, 15 == a + e + i, 15 == c + e + g} In[7]:= vars = Flatten[mat] Out[7]= {a, b, c, d, e, f, g, h, i} lin is 8 (but really 7) eqns, 9 unknowns. How to impose distinct positive integers? Hint, nonlinearly: In[8]:= cub = Total[vars^3] == Sum[n^3, {n, 9}] Out[8]= a^3 + b^3 + c^3 + d^3 + e^3 + f^3 + g^3 + h^3 + i^3 == 2025 In[9]:= quad = Total[vars^2] == Sum[n^2, {n, 9}] Out[9]= a^2 + b^2 + c^2 + d^2 + e^2 + f^2 + g^2 + h^2 + i^2 == 285 In[10]:= forth = Total[vars^4] == Sum[n^4, {n, 9}] Out[10]= a^4 + b^4 + c^4 + d^4 + e^4 + f^4 + g^4 + h^4 + i^4 == 15333 In[6]:= harm = FullSimplify[1/Total[vars^-1] == 1/Sum[n^-1, {n, 9}]] Out[6]= 1/(1/a + 1/b + 1/c + 1/d + 1/e + 1/f + 1/g + 1/h + 1/i) == 2520/7129 (Note, Solve boggles if this equation is reciprocated.) In[15]:= fac = a*b*c*d*e*f*g*h*i == 9! Out[15]= a b c d e f g h i == 362880 etc. E.g.,we get the eight sol'ns if we impose the cubic and quintic, In[15]:= tim[Solve[Join[lin, {quin, cub}]]] During evaluation of In[15]:= 0.032546 Out[15]= {{b -> 1, d -> 3, f -> 7, h -> 9, a -> 8, c -> 6, g -> 4, i -> 2, e -> 5}, {b -> 1, d -> 7, f -> 3, h -> 9, a -> 6, c -> 8, g -> 2, i -> 4, e -> 5}, {b -> 3, d -> 1, f -> 9, h -> 7, a -> 8, c -> 4, g -> 6, i -> 2, e -> 5}, {b -> 3, d -> 9, f -> 1, h -> 7, a -> 4, c -> 8, g -> 2, i -> 6, e -> 5}, {b -> 7, d -> 1, f -> 9, h -> 3, a -> 6, c -> 2, g -> 8, i -> 4, e -> 5}, {b -> 7, d -> 9, f -> 1, h -> 3, a -> 2, c -> 6, g -> 4, i -> 8, e -> 5}, {b -> 9, d -> 3, f -> 7, h -> 1, a -> 4, c -> 2, g -> 8, i -> 6, e -> 5}, {b -> 9, d -> 7, f -> 3, h -> 1, a -> 2, c -> 4, g -> 6, i -> 8, e -> 5}} tim[Solve[Join[lin, {sex, nine}]]] 3.5 sec gives the above, plus 40 solutions with ghastly surds. MinimalPolynomial boggles. Except e always ->5. tim[Solve[Join[lin, {fac, quad}]]] 0.176237 gives the msq sol'ns plus 8 with docile quartic surds. Earlier, we were puzzled that Solve[Join[lin, {quad,cub}]] remained underconstrained: In[29]:= tim[Solve[Join[lin, {cub, quad}]]] During evaluation of In[29]:= 0.020920 Out[29]= {{b -> i - Sqrt[-15 + 10 i - i^2], d -> i + Sqrt[-15 + 10 i - i^2], f -> 10 - i - Sqrt[-15 + 10 i - i^2], h -> 10 - i + Sqrt[-15 + 10 i - i^2], a -> 10 - i, c -> 5 + Sqrt[-15 + 10 i - i^2], g -> 5 - Sqrt[-15 + 10 i - i^2], e -> 5}, {b -> i + Sqrt[-15 + 10 i - i^2], d -> i - Sqrt[-15 + 10 i - i^2], f -> 10 - i + Sqrt[-15 + 10 i - i^2], h -> 10 - i - Sqrt[-15 + 10 i - i^2], a -> 10 - i, c -> 5 - Sqrt[-15 + 10 i - i^2], g -> 5 + Sqrt[-15 + 10 i - i^2], e -> 5}} (Note two sol'ns vs eight.) An "explanation": In[23]:= tim[GroebnerBasis[Subtract @@ # & /@ Join[lin, {quad}], vars]] During evaluation of In[23]:= 0.005090 Out[23]= {115 - 20 h + h^2 - 30 i + 2 h i + 2 i^2, -15 + g + h + i, -20 + f + h + 2 i, -5 + e, 10 + d - h - 2 i, 5 + c - h - i, -10 + b + h, -10 + a + i} In[24]:= tim[GroebnerBasis[Subtract @@ # & /@ Join[lin, {cub}], vars]] During evaluation of In[24]:= 0.003059 Out[24]= {115 - 20 h + h^2 - 30 i + 2 h i + 2 i^2, -15 + g + h + i, -20 + f + h + 2 i, -5 + e, 10 + d - h - 2 i, 5 + c - h - i, -10 + b + h, -10 + a + i} I.e., cub and quad are completely interchangeable in this problem. NeilB says this is because sum(n^3) = sum(n)^2. Sounds plausible, but I need clarification. For the 4x4 msqs, 4+4+2-1 eqns, 16 unknowns leaves too many nonlinear constraints for current technology, even if we add two more linears for the corners and central 2x2. --rwg