z3 is satisfiability modulo theories (SMT), docs
solve simple linear equation system like this:
import z3

solver = z3.Solver()
x = z3.Int('x') # 2
y = z3.Int('y') # 3
z = z3.Int('z') # 5

solver.add(x + y == 5)
solver.add(y + z == 8)
solver.add(z + x == 7)

if solver.check() == z3.sat:
    m = solver.model()
    print("x=" + str(m.eval(x).as_long()))
    print("y=" + str(m.eval(y).as_long()))
    print("z=" + str(m.eval(z).as_long()))
else:
    print("No solution!")
It can also solve quadratic equations. If there are several solutions, z3 will
just return one of them. If you want to know if there is just a single solution,
or if you want to find the other solutions, you can add a further constraint
saying that the solution may not be X (where X is the first solution z3 provided)
from z3 import *

x = Real('x')

solver = Solver()
solver.add(x * x - 5*x + 6 == 0)

if solver.check() == sat:
    model = solver.model()
    print("solution=" + str(model[x]))
    solver.add(x != model[x])
    if solver.check() == sat:
        model = solver.model()
        print("another solution x=" + str(model[x]))
    else:
        print("there was just one solution")
else:
    print("no solution")
A variable is called "under-constrained" if there isn't enough constraints.
If you eval an under-constrained variable it returns the "symbolic" variable
rather than concrete value(s), so you can check if a variable is under-constrained
by doing:
    if m.eval(foo, model_completion=False) == foo:
        print("variable foo is under-constrained")
if a cipher text is constructed by xoring a plaintext and a keystream, then
z3 can be used to recover the plaintext
import z3

solver = z3.Solver()

#pt = bytes.fromhex("68656c6c6f20776f726c64")
k  = bytes.fromhex("3132333435363738393031")
ct = bytes.fromhex("59575f585a1640574b5c55")
pt = [z3.BitVec(f"pt_{i}", 8) for i in range(len(ct))]
foo = z3.Int("foo")

for idx in range(0, len(ct)):
    solver.add(pt[idx] ^ k[idx] == ct[idx])

if solver.check() == z3.sat:
    m = solver.model()

    if m.eval(foo, model_completion=False) == foo:
        print("variable foo is under-constrained")

    if z3.is_bv_value(m.eval(pt[0], model_completion=False)):
        print("variable pt[0] is NOT under-constrained: " + str(m.eval(pt[0], model_completion=False).as_long()))
    else:
        print("variable pt[0] is under-constrained")

    if z3.is_bv_value(m.eval(pt[1], model_completion=False)):
        print("variable pt[1] is NOT under-constrained: " + str(m.eval(pt[1], model_completion=False).as_long()))
    else:
        print("variable pt[1] is under-constrained")


    pt_str = ''.join([chr(m.eval(pt[idx]).as_long()) for idx in range(0, len(ct))])
    print("pt=" + pt_str + "####")

else:
    print("No solution!")