If there are a finite number of solutions, you can use the disjunct of the constants (your x_i's) not equal to their assigned model values to enumerate all of them. If there are infinite solutions (which is the case if you want to prove this for all natural numbers n), you can use the same technique, but of course couldn't enumerate them all, but could use this to generate many solutions up to some bound you pick. If you want to prove this for all n > 1, you will need to use quantifiers. I've added a discussion of this below.
While you didn't quite ask this question, you should see this question/answer as well: Z3: finding all satisfying models
Here's your example doing this (z3py link here: http://rise4fun.com/Z3Py/643M ):
# Add the equations { x_1^2 == 1, x_2 == 1, ... x_n == 1 } to s and return it.
def add_constraints(s, n, model):
assert n > 1
X = IntVector('x', n)
s.add(X[0]*X[0] == 1)
for i in xrange(1, n):
s.add(X[i] == 1)
notAgain = []
i = 0
for val in model:
notAgain.append(X[i] != model[val])
i = i + 1
if len(notAgain) > 0:
s.add(Or(notAgain))
print Or(notAgain)
return s
for n in range(2,5):
s = Solver()
i = 0
add_constraints(s, n, [])
while s.check() == sat:
print s.model()
i = i + 1
add_constraints(s, n, s.model())
print i # solutions
If you want to prove there are no other solutions for any choice of n, you need to use quantifiers, since the previous approach will only work for finite n (and it gets very expensive quickly). Here is an encoding showing this proof. You could generalize this to incorporate the model generation capability in the previous part to come up with the +/- 1 solution for a more general formula. If the equation has a number of solutions independent of n (like in your example), this would allow you to prove equations have some finite number of solutions. If the number of solutions is a function of n, you'd have to figure that function out. z3py link: http://rise4fun.com/Z3Py/W9En
x = Function('x', IntSort(), IntSort())
s = Solver()
n = Int('n')
# theorem says that x(1)^2 == 1 and that x(1) != +/- 1, and forall n >= 2, x(n) == 1
# try removing the x(1) != +/- constraints
theorem = ForAll([n], And(Implies(n == 1, And(x(n) * x(n) == 1, x(n) != 1, x(n) != -1) ), Implies(n > 1, x(n) == 1)))
#s.add(Not(theorem))
s.add(theorem)
print s.check()
#print s.model() # unsat, no model available, no other solutions