Repository: sympy/sympy

Issue:
linear_eq_to_matrix: fails on simple example
SymPy 1.3 is wrong
```
>>> eqn = Eq(a*(x + 2*y) + b*y, 3)
a⋅(x + 2⋅y) + b⋅y = 3

>>> eqn_expand = eqn.expand()
a⋅x + 2⋅a⋅y + b⋅y = 3

>>> linear_eq_to_matrix([eqn], [x, y])   # WRONG
([0  b], [3])

>>> linear_eq_to_matrix([eqn_expand], [x, y])    # RIGHT
([a  2⋅a + b], [3])
```


SymPy 1.4-rc2 is a different kind of wrong:
```
>>> linear_eq_to_matrix([eqn], [x, y])      # WRONG
([a  b], [3])

>>> linear_eq_to_matrix([eqn_expand], [x, y])   # RIGHT
([a  2⋅a + b], [3])
```

@asmeurer this is sort of a regression, although it was broken in 1.3 as well.  Do you want to hold 1.4 for this or not?

Some effort went into `linear_eq_to_matrix` in this last cycle so b/c of that I'd be inclined to try to fix this; although not sure I'll get to it today...

Discussion:

Patch:
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index ad2df78a5f..db89a9fa26 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -1858,8 +1858,8 @@ def linear_coeffs(eq, *syms, **_kw):
         elif f.is_Add:
             d1 = linear_coeffs(f, *syms, **{'dict': True})
             d[0].append(m*d1.pop(0))
-            xf, vf = list(d1.items())[0]
-            d[xf].append(m*vf)
+            for xf, vf in d1.items():
+                d[xf].append(m*vf)
         else:
             break
     else:


Test patch:
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 33d7954e92..61f15fa6f9 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1060,6 +1060,11 @@ def test_linear_eq_to_matrix():
         Matrix([[1]]), Matrix([[-1]]))
 
 
+def test_issue_16577():
+    assert linear_eq_to_matrix(Eq(a*(2*x + 3*y) + 4*y, 5), x, y) == (
+        Matrix([[2*a, 3*a + 4]]), Matrix([[5]]))
+
+
 def test_linsolve():
     x, y, z, u, v, w = symbols("x, y, z, u, v, w")
     x1, x2, x3, x4 = symbols('x1, x2, x3, x4')
@@ -2025,3 +2030,4 @@ def test_linear_coeffs():
         linear_coeffs(x + 2*x**2 + x**3, x, x**2))
     raises(ValueError, lambda:
         linear_coeffs(1/x*(x - 1) + 1/x, x))
+    assert linear_coeffs(a*(x + y), x, y) == [a, a, 0]