Presumably if you have 20 million equations, they came from a program that's has fewer than 20 million moving parts, like if they came from A x = b where the matrix A has 20M entries. The gist is either exploit structure to make a massive number of small equations or keep the symbols in their "natural" form instead of reducing to scalars, and work with more advanced CAS functionality (like, you might have to learn about noncommutative variations on groebner bases). But also, yes sympy is ultra slow with some things.