1. If your system is sparse, meaning most of the coefficients are zeroes, you can exploit this with any decomposition (factorization) or elimination based method that works on sparse structures.
  2. If your matrix is diagonally dominant, meaning that the largest elements in every equation tend to land on the diagonal, then iterative methods such as Jacobi or Gauss-Seidel should work best for you.
  3. And if your system is very small, like 2 to 4 equations, you're better off without any solver at all. Just solve it symbolically with a SymPy and copy-paste the symbolic solution back to your code.
←
→