A min/max finder based on Lagrange Multipliers with one constraint

lagrange-solver-one-constraint.py

from sympy import *
from sympy.abc import x, y, z, l

f = x**2 + y**2 + z**2
g = x**4 + y**4 + z**4
k = 9

grad = lambda fn: Matrix([fn.diff(x), fn.diff(y), fn.diff(z)]).transpose()

# grad(f) = lambda * grad(g)
first = [*(grad(f) - grad(g) * l)]
second = g - k

solns = solve([*first, second], (x, y, z, l))
final_solns = []
for soln in solns:
    # filter out solutions with imaginary parts
    if not (soln[0].is_imaginary or soln[1].is_imaginary or soln[2].is_imaginary or soln[3].is_imaginary):
        final_solns += [soln]

vals = [f.subs({'x': x[0], 'y': x[1], 'z': x[2] }) for x in final_solns]

print('max value: ', max(vals))
print('min value: ', min(vals))