Check out SymPy:

    >>> x, a, b, c = symbols('x a b c')
    >>> solve(a*x**2+b*x+c, x)     # solve(expr_eqal_to_zero, var)
    [ (-b + sqrt(-4*a*c + b**2))/(2*a),
      (-b - sqrt(-4*a*c + b**2))/(2*a)  ]

More info: https://minireference.com/static/tutorials/sympy_tutorial.pd...