Is there software that helps list small solutions of the Diophantine equation $$ x_0^2=1+x_1^2+x_2^2+\cdots+ x_n^2 $$ where "small" is negotiable, but e.g. we could fix $x_0$ and and ask for the list of all solutions $x_1, \dots, x_n\in\mathbb Z$?

2$\begingroup$ This is $$(x_0x_1)(x_0+x_1)=x_0^2x_1^2=1+x_2^2+\cdots +x_n^2,$$ so sort of reduces to factorization. $\endgroup$ – Robin Chapman Jun 10 '10 at 20:17

$\begingroup$ @Will, I wish to distinguish $x_j$ and $x_j$. Please do not bother to write anything in C++; I just want to know whether there is anything already available as part of Mathematica ot other similar package. @Robin, thanks, reducing factorization is helpful idea. $\endgroup$ – Igor Belegradek Jun 10 '10 at 21:01

2$\begingroup$ For $n=2$ you can try alpertron.com.ar/QUAD.HTM for various values of $x_0$. $\endgroup$ – Eric Rowell Jun 10 '10 at 21:12

2$\begingroup$ SquaresR[d, n] gives the number of ways r_d (n) to represent the integer n as a sum of d squares. PowersRepresentations[n, k, p] gives the distinct representations of the integer n as a sum of k nonnegative p\^th integer powers. EllipticTheta[a,u,q] gives the theta function Subscript[[CurlyTheta], a](u,q) (a=1,[Ellipsis],4). $\endgroup$ – Will Jagy Jun 10 '10 at 21:19

$\begingroup$ note the above are Mathematica commands $\endgroup$ – Will Jagy Jun 10 '10 at 21:52
Have you had a look on this tutorial: http://reference.wolfram.com/mathematica/tutorial/DiophantineReduce.html ? To take different values n < n_max into account, a simple loop could work.
This function in Mathematica find them all:
f[n_] := Reduce[ Total[Table[x[i]^2, {i, 1, n}]] + 1 == x[0]^2 &&
(Table[x[i], {i, 1, 3}] /. List > LessEqual),
Table[x[i], {i, 0, n}], Integers]
Invoke with:
f[1], f[3] ... etc