I am trying to solve the equation, where lambda is a (possibly complex) scalar and A, B, and C are N x N matrices. The goal is to compute the roots of the characteristic polynomial in lambda. I have successfully done this for a low dimensional system using sympy (N = 4). However, the calculation hangs for N as small as 10. With A, B, and C initialized as numpy arrays, I use the following code: c = sympy.Symbol('c') # Define constant mat = c*c*sympy.Matrix(A) + c*sympy.Matrix(B) + sympy.Matrix(C) # Assemble matrix det = mat.det(method='berkowitz') # Compute determinant soln = sympy.solve(det,c,dict=True) # Find roots The calculation hangs at the third line. The end goal is to do this computation for N ~ 1000, so I am wondering how I can speed up the calculation. I tried computing |A + B + C| for N ~ 500 and this took less than a second, so it's something about the symbolic calculation that increases the time. Additionally, similar computations in Mathematica seem to have no issues. Does the mat.det() function treat all matrix entries as symbols? All of the elements of A, B, and C are floating point, so I am also wondering if there is a way to treat all coefficients except lambda numerically. Continue reading...