Thanks Robert, I enjoyed reading this account.
Two notes, one mathematical and one Maple-related. First the mathematical one: I think you would expect 3 degrees of freedom, right? There are 10 equations, but there is a dependency between them: all four row-sums must needs sum to the same as all four column-sums. So there are really only 9 independent equations.
The Maple point is that the combinatorial question you raise is solved by the Iterator package, in particular the BoundedComposition command. You supply it with a k-tuple of bounds (b1, ..., bk), and a sum s, all nonnegative integers. It then finds all k-tuples (a1, ..., ak) of nonnegative integers with ai <= bi that sum to s. This raises one problem: we want integers from 1-9, and BoundedComposition will also give us zeroes -- but that's easily dealt with: we just use integers from 0-8 instead, adjust the requested sum (by subtracting k=4), and presto. For example:
k := 4;
B := BoundedComposition([8 $ k], 25-k):
Number(B); # result: 324
tells us that there are 324 such tuples with sum 25, and will print all of them. To do a computation with them, you would use:
for b in BoundedComposition([8 $ k], 25 - k) do
tup := b +~ 1;
# Now use the tuple, tup, which will contain suitable integers 1-9. (It's an Array.)
Of course, if one number is given already, you could generate 3-tuples.