Repository: sympy/sympy

Issue:
Wrong linear_eq_to_matrix if equations not expanded
I noticed that `linear_eq_to_matrix` doesn't yield the correct matrix if the linear equations are not fully expanded. For example
```python
>>> x1, x2, mu = symbols('x1, x2, mu')
>>> equations = Matrix([[mu * (x1 + x2) + mu],
>>>                     [5 * x2 - 1]])
>>> A, b = linear_eq_to_matrix(equations, symbols)
>>> A
```
gives the wrong result `Matrix([[0, 0], [0, 5]])`. This can be fixed by adding `equations = equations.expand()` prior to the `linear_eq_to_matrix` call. Would it make sense to include an `expand` call in the `linear_eq_to_matrix` method itself? Willing to submit the pull request.

Discussion:

Patch:
diff --git a/sympy/matrices/common.py b/sympy/matrices/common.py
index 98704cc91f..af3c6fb25f 100644
--- a/sympy/matrices/common.py
+++ b/sympy/matrices/common.py
@@ -437,6 +437,7 @@ def row_insert(self, pos, other):
         row
         col_insert
         """
+        from sympy.matrices import MutableMatrix
         # Allows you to build a matrix even if it is null matrix
         if not self:
             return self._new(other)
diff --git a/sympy/physics/continuum_mechanics/beam.py b/sympy/physics/continuum_mechanics/beam.py
index 8d11dbe859..387c105d1f 100644
--- a/sympy/physics/continuum_mechanics/beam.py
+++ b/sympy/physics/continuum_mechanics/beam.py
@@ -1490,13 +1490,13 @@ class Beam3D(Beam):
     >>> b.bc_deflection = [(0, [0, 0, 0]), (l, [0, 0, 0])]
     >>> b.solve_slope_deflection()
     >>> b.slope()
-    [0, 0, l*x*(-l*q + 3*l*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I)) + 3*m)/(6*E*I)
-    + q*x**3/(6*E*I) + x**2*(-l*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I))
+    [0, 0, l*x*(-l*q + 3*l*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I)) + 3*m)/(6*E*I)
+    + q*x**3/(6*E*I) + x**2*(-l*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*(A*G*l**2 + 12*E*I))
     - m)/(2*E*I)]
     >>> b.deflection()
-    [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I))
-    + l*m*x**2/(4*E*I) - l*x**3*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I)
-    + q*x**4/(24*E*I) + l*x*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I)) - q*x**2/(2*A*G), 0]
+    [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I))
+    + l*m*x**2/(4*E*I) - l*x**3*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I)
+    + q*x**4/(24*E*I) + l*x*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I)) - q*x**2/(2*A*G), 0]
 
     References
     ==========
diff --git a/sympy/release.py b/sympy/release.py
index 57dc05521c..80b14a44ed 100644
--- a/sympy/release.py
+++ b/sympy/release.py
@@ -1,1 +1,1 @@
-__version__ = "1.4.dev"
+__version__ = "1.2.1.dev"
diff --git a/sympy/series/sequences.py b/sympy/series/sequences.py
index 1959e7ac62..08f36df061 100644
--- a/sympy/series/sequences.py
+++ b/sympy/series/sequences.py
@@ -3,8 +3,7 @@
 from sympy.core.basic import Basic
 from sympy.core.mul import Mul
 from sympy.core.singleton import S, Singleton
-from sympy.core.symbol import Dummy, Symbol, Wild
-from sympy.core.function import UndefinedFunction
+from sympy.core.symbol import Dummy, Symbol
 from sympy.core.compatibility import (range, integer_types, with_metaclass,
                                       is_sequence, iterable, ordered)
 from sympy.core.decorators import call_highest_priority
@@ -12,8 +11,6 @@
 from sympy.core.sympify import sympify
 from sympy.core.containers import Tuple
 from sympy.core.evaluate import global_evaluate
-from sympy.core.numbers import Integer
-from sympy.core.relational import Eq
 from sympy.polys import lcm, factor
 from sympy.sets.sets import Interval, Intersection
 from sympy.utilities.iterables import flatten
@@ -714,173 +711,6 @@ def coeff_mul(self, coeff):
         return SeqFormula(formula, self.args[1])
 
 
-class RecursiveSeq(SeqBase):
-    """A finite degree recursive sequence.
-
-    That is, a sequence a(n) that depends on a fixed, finite number of its
-    previous values. The general form is
-
-        a(n) = f(a(n - 1), a(n - 2), ..., a(n - d))
-
-    for some fixed, positive integer d, where f is some function defined by a
-    SymPy expression.
-
-    Parameters
-    ==========
-
-    recurrence : SymPy expression defining recurrence
-        This is *not* an equality, only the expression that the nth term is
-        equal to. For example, if :code:`a(n) = f(a(n - 1), ..., a(n - d))`,
-        then the expression should be :code:`f(a(n - 1), ..., a(n - d))`.
-
-    y : function
-        The name of the recursively defined sequence without argument, e.g.,
-        :code:`y` if the recurrence function is :code:`y(n)`.
-
-    n : symbolic argument
-        The name of the variable that the recurrence is in, e.g., :code:`n` if
-        the recurrence function is :code:`y(n)`.
-
-    initial : iterable with length equal to the degree of the recurrence
-        The initial values of the recurrence.
-
-    start : start value of sequence (inclusive)
-
-    Examples
-    ========
-
-    >>> from sympy import Function, symbols
-    >>> from sympy.series.sequences import RecursiveSeq
-    >>> y = Function("y")
-    >>> n = symbols("n")
-    >>> fib = RecursiveSeq(y(n - 1) + y(n - 2), y, n, [0, 1])
-
-    >>> fib.coeff(3) # Value at a particular point
-    2
-
-    >>> fib[:6] # supports slicing
-    [0, 1, 1, 2, 3, 5]
-
-    >>> fib.recurrence # inspect recurrence
-    Eq(y(n), y(n - 2) + y(n - 1))
-
-    >>> fib.degree # automatically determine degree
-    2
-
-    >>> for x in zip(range(10), fib): # supports iteration
-    ...     print(x)
-    (0, 0)
-    (1, 1)
-    (2, 1)
-    (3, 2)
-    (4, 3)
-    (5, 5)
-    (6, 8)
-    (7, 13)
-    (8, 21)
-    (9, 34)
-
-    See Also
-    ========
-
-    sympy.series.sequences.SeqFormula
-
-    """
-
-    def __new__(cls, recurrence, y, n, initial=None, start=0):
-        if not isinstance(y, UndefinedFunction):
-            raise TypeError("recurrence sequence must be an undefined function"
-                            ", found `{}`".format(y))
-
-        if not isinstance(n, Basic) or not n.is_symbol:
-            raise TypeError("recurrence variable must be a symbol"
-                            ", found `{}`".format(n))
-
-        k = Wild("k", exclude=(n,))
-        degree = 0
-
-        # Find all applications of y in the recurrence and check that:
-        #   1. The function y is only being used with a single argument; and
-        #   2. All arguments are n + k for constant negative integers k.
-
-        prev_ys = recurrence.find(y)
-        for prev_y in prev_ys:
-            if len(prev_y.args) != 1:
-                raise TypeError("Recurrence should be in a single variable")
-
-            shift = prev_y.args[0].match(n + k)[k]
-            if not (shift.is_constant() and shift.is_integer and shift < 0):
-                raise TypeError("Recurrence should have constant,"
-                                " negative, integer shifts"
-                                " (found {})".format(prev_y))
-
-            if -shift > degree:
-                degree = -shift
-
-        if not initial:
-            initial = [Dummy("c_{}".format(k)) for k in range(degree)]
-
-        if len(initial) != degree:
-            raise ValueError("Number of initial terms must equal degree")
-
-        degree = Integer(degree)
-        start = sympify(start)
-
-        initial = Tuple(*(sympify(x) for x in initial))
-
-        seq = Basic.__new__(cls, recurrence, y(n), initial, start)
-
-        seq.cache = {y(start + k): init for k, init in enumerate(initial)}
-        seq._start = start
-        seq.degree = degree
-        seq.y = y
-        seq.n = n
-        seq._recurrence = recurrence
-
-        return seq
-
-    @property
-    def start(self):
-        """The starting point of the sequence. This point is included"""
-        return self._start
-
-    @property
-    def stop(self):
-        """The ending point of the sequence. (oo)"""
-        return S.Infinity
-
-    @property
-    def interval(self):
-        """Interval on which sequence is defined."""
-        return (self._start, S.Infinity)
-
-    def _eval_coeff(self, index):
-        if index - self._start < len(self.cache):
-            return self.cache[self.y(index)]
-
-        for current in range(len(self.cache), index + 1):
-            # Use xreplace over subs for performance.
-            # See issue #10697.
-            seq_index = self._start + current
-            current_recurrence = self._recurrence.xreplace({self.n: seq_index})
-            new_term = current_recurrence.xreplace(self.cache)
-
-            self.cache[self.y(seq_index)] = new_term
-
-        return self.cache[self.y(self._start + current)]
-
-    def __iter__(self):
-        index = self._start
-        while True:
-            yield self._eval_coeff(index)
-            index += 1
-
-    @property
-    def recurrence(self):
-        """Equation defining recurrence."""
-        return Eq(self.y(self.n), self._recurrence)
-
-
 def sequence(seq, limits=None):
     """Returns appropriate sequence object.
 
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index c0a243b3d9..281d06b8ce 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -19,7 +19,7 @@
 from sympy.core.facts import InconsistentAssumptions
 from sympy.core.numbers import I, Number, Rational, oo
 from sympy.core.function import (Lambda, expand_complex, AppliedUndef,
-                                expand_log)
+                                expand_log, _mexpand)
 from sympy.core.relational import Eq, Ne
 from sympy.core.symbol import Symbol
 from sympy.simplify.simplify import simplify, fraction, trigsimp
@@ -34,7 +34,7 @@
 from sympy.sets import (FiniteSet, EmptySet, imageset, Interval, Intersection,
                         Union, ConditionSet, ImageSet, Complement, Contains)
 from sympy.sets.sets import Set
-from sympy.matrices import Matrix
+from sympy.matrices import Matrix, MatrixBase
 from sympy.polys import (roots, Poly, degree, together, PolynomialError,
                          RootOf, factor)
 from sympy.solvers.solvers import (checksol, denoms, unrad,
@@ -42,11 +42,12 @@
 from sympy.solvers.polysys import solve_poly_system
 from sympy.solvers.inequalities import solve_univariate_inequality
 from sympy.utilities import filldedent
-from sympy.utilities.iterables import numbered_symbols
+from sympy.utilities.iterables import numbered_symbols, has_dups
 from sympy.calculus.util import periodicity, continuous_domain
 from sympy.core.compatibility import ordered, default_sort_key, is_sequence
 
 from types import GeneratorType
+from collections import defaultdict
 
 
 def _masked(f, *atoms):
@@ -1806,13 +1807,73 @@ def solvify(f, symbol, domain):
 ###############################################################################
 
 
+def linear_coeffs(eq, *syms):
+    """Return a dictionary with symbols in `syms` keyed to the
+    coefficient of the corresponding symbol in `eq`. The additive
+    constant is returned with key 0.
+
+    Examples
+    ========
+
+    >>> from sympy.solvers.solveset import linear_coeffs
+    >>> from sympy.abc import x, y, z
+
+
+    >>> linear_coeffs(2*x - 1, x, y)
+    {0: -1, x: 2}
+
+    It is not necessary to expand the expression:
+
+    >>> linear_coeffs(x + y*(z*(3*x + 2) + 3), x)
+    {0: y*(2*z + 3), x: 3*y*z + 1}
+
+    But if there are nonlinear or cross terms -- even if they would
+    cancel after simplification -- an error is raised so the situation
+    does not pass silently past the caller's attention:
+
+    >>> eq = 1/x*(x - 1) + 1/x
+    >>> linear_coeffs(eq.expand(), x)
+    {0: 1}
+    >>> linear_coeffs(eq, x)
+    Traceback (most recent call last):
+    ...
+    ValueError: nonlinear term encountered: 1/x
+
+    >>> linear_coeffs(x*(y + 1) - x*y, x, y)
+    Traceback (most recent call last):
+    ...
+    ValueError: nonlinear term encountered: x*(y + 1)
+    """
+    d = defaultdict(list)
+    c, terms = eq.as_coeff_add(*syms)
+    d[0].append(c)
+    for t in terms:
+        m, f = t.as_coeff_mul(*syms)
+        if len(f) != 1:
+            break
+        f = f[0]
+        if f in syms:
+            d[f].append(m)
+        elif f.is_Add:
+            d1 = linear_coeffs(f, *syms)
+            d[0].append(m*d1.pop(0))
+            xf, vf = list(d1.items())[0]
+            d[xf].append(m*vf)
+        else:
+            break
+    else:
+        for k, v in d.items():
+            d[k] = Add(*v)
+        return defaultdict(lambda: S.Zero, d)  # default now zero
+    raise ValueError('nonlinear term encountered: %s' % t)
+
+
 def linear_eq_to_matrix(equations, *symbols):
     r"""
     Converts a given System of Equations into Matrix form.
     Here `equations` must be a linear system of equations in
-    `symbols`. The order of symbols in input `symbols` will
-    determine the order of coefficients in the returned
-    Matrix.
+    `symbols`. Element M[i, j] corresponds to the coefficient
+    of the jth symbol in the ith equation.
 
     The Matrix form corresponds to the augmented matrix form.
     For example:
@@ -1829,84 +1890,102 @@ def linear_eq_to_matrix(equations, *symbols):
      A = [ 3  1  1 ]   b  =   [-6 ]
          [ 2  4  9 ]          [ 2 ]
 
+    The only simplification performed is to convert
+    `Eq(a, b) -> a - b`.
+
+    Raises
+    ======
+
+    ValueError
+        The equations contain a nonlinear term.
+        The symbols are not given or are not unique.
+
     Examples
     ========
 
     >>> from sympy import linear_eq_to_matrix, symbols
-    >>> x, y, z = symbols('x, y, z')
-    >>> eqns = [x + 2*y + 3*z - 1, 3*x + y + z + 6, 2*x + 4*y + 9*z - 2]
-    >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
-    >>> A
-    Matrix([
-    [1, 2, 3],
-    [3, 1, 1],
-    [2, 4, 9]])
-    >>> b
-    Matrix([
-    [ 1],
-    [-6],
-    [ 2]])
-    >>> eqns = [x + z - 1, y + z, x - y]
+    >>> c, x, y, z = symbols('c, x, y, z')
+
+    The coefficients (numerical or symbolic) of the symbols will
+    be returned as matrices:
+
+    >>> eqns = [c*x + z - 1 - c, y + z, x - y]
     >>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
     >>> A
     Matrix([
-    [1,  0, 1],
+    [c,  0, 1],
     [0,  1, 1],
     [1, -1, 0]])
     >>> b
     Matrix([
-    [1],
-    [0],
-    [0]])
-
-    * Symbolic coefficients are also supported
-
-    >>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
-    >>> eqns = [a*x + b*y - c, d*x + e*y - f]
-    >>> A, B = linear_eq_to_matrix(eqns, x, y)
-    >>> A
-    Matrix([
-    [a, b],
-    [d, e]])
-    >>> B
-    Matrix([
-    [c],
-    [f]])
-
+    [c + 1],
+    [    0],
+    [    0]])
+
+    This routine does not simplify expressions and will raise an error
+    if nonlinearity is encountered:
+
+    >>> eqns = [
+    ...     (x**2 - 3*x)/(x - 3) - 3,
+    ...     y**2 - 3*y - y*(y - 4) + x - 4]
+    >>> linear_eq_to_matrix(eqns, [x, y])
+    Traceback (most recent call last):
+    ...
+    ValueError:
+    The term (x**2 - 3*x)/(x - 3) is nonlinear in {x, y}
+
+    Simplifying these equations will discard the removable singularity
+    in the first, reveal the linear structure of the second:
+
+    >>> [e.simplify() for e in eqns]
+    [x - 3, x + y - 4]
+
+    Any such simplification needed to eliminate nonlinear terms must
+    be done before calling this routine.
     """
-
     if not symbols:
-        raise ValueError('Symbols must be given, for which coefficients \
-                         are to be found.')
+        raise ValueError(filldedent('''
+            Symbols must be given, for which coefficients
+            are to be found.
+            '''))
 
     if hasattr(symbols[0], '__iter__'):
         symbols = symbols[0]
 
-    M = Matrix([symbols])
-    # initialize Matrix with symbols + 1 columns
-    M = M.col_insert(len(symbols), Matrix([1]))
-    row_no = 1
+    for i in symbols:
+        if not isinstance(i, Symbol):
+            raise ValueError(filldedent('''
+            Expecting a Symbol but got %s
+            ''' % i))
+
+    if has_dups(symbols):
+        raise ValueError('Symbols must be unique')
+
+    equations = sympify(equations)
+    if isinstance(equations, MatrixBase):
+        equations = list(equations)
+    elif isinstance(equations, Expr):
+        equations = [equations]
+    elif not is_sequence(equations):
+        raise ValueError(filldedent('''
+            Equation(s) must be given as a sequence, Expr,
+            Eq or Matrix.
+            '''))
 
-    for equation in equations:
-        f = sympify(equation)
+    A, b = [], []
+    for i, f in enumerate(equations):
         if isinstance(f, Equality):
-            f = f.lhs - f.rhs
-
-        # Extract coeff of symbols
-        coeff_list = []
-        for symbol in symbols:
-            coeff_list.append(f.coeff(symbol))
-
-        # append constant term (term free from symbols)
-        coeff_list.append(-f.as_coeff_add(*symbols)[0])
-
-        # insert equations coeff's into rows
-        M = M.row_insert(row_no, Matrix([coeff_list]))
-        row_no += 1
-
-    # delete the initialized (Ist) trivial row
-    M.row_del(0)
-    A, b = M[:, :-1], M[:, -1:]
+            # non-cancelling method of creating expr; we
+            # don't care about the order of the args or
+            # else we would use _unevaluated_Add
+            f = Add._from_args(
+                Add.make_args(f.lhs) +
+                Add.make_args(-f.rhs))
+        co = linear_coeffs(f, *symbols)
+        coeff_list = [co[i] for i in symbols]
+        b.append(-co[0])
+        A.append(coeff_list)
+    A, b = map(Matrix, (A, b))
     return A, b
 
 
@@ -2053,6 +2132,16 @@ def linsolve(system, *symbols):
     >>> linsolve([], x)
     EmptySet()
 
+    * An error is raised if, after expansion, any nonlinearity
+      is detected:
+
+    >>> linsolve([x*(1/x - 1), (y - 1)**2 - y**2 + 1], x, y)
+    {(1, 1)}
+    >>> linsolve([x**2 - 1], x)
+    Traceback (most recent call last):
+    ...
+    ValueError:
+    The term x**2 is nonlinear in {x}
     """
     if not system:
         return S.EmptySet
@@ -2082,19 +2171,9 @@ def linsolve(system, *symbols):
                     symbols for which a solution is being sought must
                     be given as a sequence, too.
                 '''))
-            system = list(system)
-            for i, eq in enumerate(system):
-                try:
-                    # since we are checking it, we might as well take the
-                    # expanded expr that it will give
-                    system[i] = Poly(eq, symbols).as_expr()
-                    if any (degree(eq, sym) > 1 for sym in symbols):
-                        raise PolynomialError
-                except PolynomialError:
-                    raise ValueError(filldedent('''
-                        %s contains non-linear terms in the
-                        variables to be evaluated
-                    ''') % eq)
+            system = [
+                _mexpand(i.lhs - i.rhs if isinstance(i, Eq) else i,
+                recursive=True) for i in system]
             system, symbols, swap = recast_to_symbols(system, symbols)
             A, b = linear_eq_to_matrix(system, symbols)
             syms_needed_msg = 'free symbols in the equations provided'


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index f0a801914c..538682f032 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -10,7 +10,7 @@
 import io
 
 from sympy import (Basic, S, symbols, sqrt, sin, oo, Interval, exp, Lambda, pi,
-                   Eq, log, Function)
+                   Eq, log)
 
 from sympy.core.compatibility import range
 from sympy.utilities.pytest import XFAIL, SKIP
@@ -3715,14 +3715,6 @@ def test_sympy__series__sequences__SeqFormula():
     assert _test_args(SeqFormula(x**2, (0, 10)))
 
 
-def test_sympy__series__sequences__RecursiveSeq():
-    from sympy.series.sequences import RecursiveSeq
-    y = Function("y")
-    n = symbols("n")
-    assert _test_args(RecursiveSeq(y(n - 1) + y(n - 2), y, n, (0, 1)))
-    assert _test_args(RecursiveSeq(y(n - 1) + y(n - 2), y, n))
-
-
 def test_sympy__series__sequences__SeqExprOp():
     from sympy.series.sequences import SeqExprOp, sequence
     s1 = sequence((1, 2, 3))
diff --git a/sympy/parsing/tests/test_autolev.py b/sympy/parsing/tests/test_autolev.py
index 9f8e7eb7b4..515ee5ddba 100644
--- a/sympy/parsing/tests/test_autolev.py
+++ b/sympy/parsing/tests/test_autolev.py
@@ -1,25 +1,20 @@
-import os
-
-from sympy import sin, cos
 from sympy.external import import_module
-from sympy.utilities.pytest import skip
-from sympy.parsing.autolev import parse_autolev
+import os
+import sys
+import itertools
 
 antlr4 = import_module("antlr4")
 
 if not antlr4:
     disabled = True
 
-FILE_DIR = os.path.dirname(
-    os.path.dirname(os.path.abspath(os.path.realpath(__file__))))
-
 
 def _test_examples(in_filename, out_filename, test_name=""):
+    from sympy.parsing.autolev import parse_autolev
 
-    in_file_path = os.path.join(FILE_DIR, 'autolev', 'test-examples',
-                                in_filename)
-    correct_file_path = os.path.join(FILE_DIR, 'autolev', 'test-examples',
-                                     out_filename)
+    dir_path = os.path.dirname(os.path.dirname(os.path.abspath(os.path.realpath(__file__))))
+    in_file_path = os.path.join(dir_path, 'autolev', 'test-examples', in_filename)
+    correct_file_path = os.path.join(dir_path, 'autolev', 'test-examples', out_filename)
     with open(in_file_path) as f:
         generated_code = parse_autolev(f, include_numeric=True)
 
@@ -28,19 +23,15 @@ def _test_examples(in_filename, out_filename, test_name=""):
             if line1.startswith("#"):
                 break
             try:
-                line2 = generated_code.split('\n')[idx]
-                assert line1.rstrip() == line2.rstrip()
+               line2 = generated_code.split('\n')[idx]
+               assert line1.rstrip() == line2.rstrip()
             except Exception:
-                msg = 'mismatch in ' + test_name + ' in line no: {0}'
-                raise AssertionError(msg.format(idx+1))
+                raise AssertionError('mismatch in ' + test_name + ' in line no: {0}'.format(idx+1))
 
 
 def test_rule_tests():
-
-    l = ["ruletest1", "ruletest2", "ruletest3", "ruletest4", "ruletest5",
-         "ruletest6", "ruletest7", "ruletest8", "ruletest9", "ruletest10",
-         "ruletest11", "ruletest12"]
-
+    l = ["ruletest1", "ruletest2", "ruletest3", "ruletest4", "ruletest5", "ruletest6",\
+         "ruletest7", "ruletest8", "ruletest9", "ruletest10", "ruletest11", "ruletest12"]
     for i in l:
         in_filepath = i + ".al"
         out_filepath = i + ".py"
@@ -48,10 +39,7 @@ def test_rule_tests():
 
 
 def test_pydy_examples():
-
-    l = ["mass_spring_damper", "chaos_pendulum", "double_pendulum",
-         "non_min_pendulum"]
-
+    l = ["mass_spring_damper", "chaos_pendulum", "double_pendulum", "non_min_pendulum"]
     for i in l:
         in_filepath = os.path.join("pydy-example-repo", i + ".al")
         out_filepath = os.path.join("pydy-example-repo", i + ".py")
@@ -59,13 +47,9 @@ def test_pydy_examples():
 
 
 def test_autolev_tutorial():
-
-    dir_path = os.path.join(FILE_DIR, 'autolev', 'test-examples',
-                            'autolev-tutorial')
-
+    dir_path = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(os.path.realpath(__file__)))), 'autolev', 'test-examples', 'autolev-tutorial')
     if(os.path.isdir(dir_path)):
-        l = ["tutor1", "tutor2", "tutor3", "tutor4", "tutor5", "tutor6",
-             "tutor7"]
+        l = ["tutor1", "tutor2", "tutor3", "tutor4", "tutor5", "tutor6", "tutor7"]
         for i in l:
             in_filepath = os.path.join("autolev-tutorial", i + ".al")
             out_filepath = os.path.join("autolev-tutorial", i + ".py")
@@ -73,106 +57,15 @@ def test_autolev_tutorial():
 
 
 def test_dynamics_online():
-
-    dir_path = os.path.join(FILE_DIR, 'autolev', 'test-examples',
-                            'dynamics-online')
-
+    dir_path = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(os.path.realpath(__file__)))), 'autolev', 'test-examples', 'dynamics-online')
     if(os.path.isdir(dir_path)):
         ch1 = ["1-4", "1-5", "1-6", "1-7", "1-8", "1-9_1", "1-9_2", "1-9_3"]
-        ch2 = ["2-1", "2-2", "2-3", "2-4", "2-5", "2-6", "2-7", "2-8", "2-9",
-               "circular"]
-        ch3 = ["3-1_1", "3-1_2", "3-2_1", "3-2_2", "3-2_3", "3-2_4", "3-2_5",
-               "3-3"]
+        ch2 = ["2-1", "2-2", "2-3", "2-4", "2-5", "2-6", "2-7", "2-8", "2-9", "circular"]
+        ch3 = ["3-1_1", "3-1_2", "3-2_1", "3-2_2", "3-2_3", "3-2_4", "3-2_5", "3-3"]
         ch4 = ["4-1_1", "4-2_1", "4-4_1", "4-4_2", "4-5_1", "4-5_2"]
         chapters = [(ch1, "ch1"), (ch2, "ch2"), (ch3, "ch3"), (ch4, "ch4")]
         for ch, name in chapters:
             for i in ch:
-                in_filepath = os.path.join("dynamics-online", name, i + ".al")
-                out_filepath = os.path.join("dynamics-online", name, i + ".py")
+                in_filepath = os.path.join("dynamics-online", name , i + ".al")
+                out_filepath = os.path.join("dynamics-online", name , i + ".py")
                 _test_examples(in_filepath, out_filepath, i)
-
-
-def test_output_01():
-    """Autolev example calculates the position, velocity, and accleration of a
-    point and expresses in a single reference frame::
-
-          (1) FRAMES C,D,F
-          (2) VARIABLES FD'',DC''
-          (3) CONSTANTS R,L
-          (4) POINTS O,E
-          (5) SIMPROT(F,D,1,FD)
-       -> (6) F_D = [1, 0, 0; 0, COS(FD), -SIN(FD); 0, SIN(FD), COS(FD)]
-          (7) SIMPROT(D,C,2,DC)
-       -> (8) D_C = [COS(DC), 0, SIN(DC); 0, 1, 0; -SIN(DC), 0, COS(DC)]
-          (9) W_C_F> = EXPRESS(W_C_F>, F)
-       -> (10) W_C_F> = FD'*F1> + COS(FD)*DC'*F2> + SIN(FD)*DC'*F3>
-          (11) P_O_E>=R*D2>-L*C1>
-          (12) P_O_E>=EXPRESS(P_O_E>, D)
-       -> (13) P_O_E> = -L*COS(DC)*D1> + R*D2> + L*SIN(DC)*D3>
-          (14) V_E_F>=EXPRESS(DT(P_O_E>,F),D)
-       -> (15) V_E_F> = L*SIN(DC)*DC'*D1> - L*SIN(DC)*FD'*D2> + (R*FD'+L*COS(DC)*DC')*D3>
-          (16) A_E_F>=EXPRESS(DT(V_E_F>,F),D)
-       -> (17) A_E_F> = L*(COS(DC)*DC'^2+SIN(DC)*DC'')*D1> + (-R*FD'^2-2*L*COS(DC)*DC'*FD'-L*SIN(DC)*FD'')*D2> + (R*FD''+L*COS(DC)*DC''-L*SIN(DC)*DC'^2-L*SIN(DC)*FD'^2)*D3>
-
-    """
-
-    if not antlr4:
-        skip('Test skipped: antlr4 is not installed.')
-
-    autolev_input = """\
-FRAMES C,D,F
-VARIABLES FD'',DC''
-CONSTANTS R,L
-POINTS O,E
-SIMPROT(F,D,1,FD)
-SIMPROT(D,C,2,DC)
-W_C_F>=EXPRESS(W_C_F>,F)
-P_O_E>=R*D2>-L*C1>
-P_O_E>=EXPRESS(P_O_E>,D)
-V_E_F>=EXPRESS(DT(P_O_E>,F),D)
-A_E_F>=EXPRESS(DT(V_E_F>,F),D)\
-"""
-
-    sympy_input = parse_autolev(autolev_input)
-
-    g = {}
-    l = {}
-    exec(sympy_input, g, l)
-
-    w_c_f = l['frame_c'].ang_vel_in(l['frame_f'])
-    # P_O_E> means "the position of point E wrt to point O"
-    p_o_e = l['point_e'].pos_from(l['point_o'])
-    v_e_f = l['point_e'].vel(l['frame_f'])
-    a_e_f = l['point_e'].acc(l['frame_f'])
-
-    # NOTE : The Autolev outputs above were manually transformed into
-    # equivalent SymPy physics vector expressions. Would be nice to automate
-    # this transformation.
-    expected_w_c_f = (l['fd'].diff()*l['frame_f'].x +
-                      cos(l['fd'])*l['dc'].diff()*l['frame_f'].y +
-                      sin(l['fd'])*l['dc'].diff()*l['frame_f'].z)
-
-    assert (w_c_f - expected_w_c_f).simplify() == 0
-
-    expected_p_o_e = (-l['l']*cos(l['dc'])*l['frame_d'].x +
-                      l['r']*l['frame_d'].y +
-                      l['l']*sin(l['dc'])*l['frame_d'].z)
-
-    assert (p_o_e - expected_p_o_e).simplify() == 0
-
-    expected_v_e_f = (l['l']*sin(l['dc'])*l['dc'].diff()*l['frame_d'].x -
-                      l['l']*sin(l['dc'])*l['fd'].diff()*l['frame_d'].y +
-                      (l['r']*l['fd'].diff() +
-                       l['l']*cos(l['dc'])*l['dc'].diff())*l['frame_d'].z)
-    assert (v_e_f - expected_v_e_f).simplify() == 0
-
-    expected_a_e_f = (l['l']*(cos(l['dc'])*l['dc'].diff()**2 +
-                              sin(l['dc'])*l['dc'].diff().diff())*l['frame_d'].x +
-                      (-l['r']*l['fd'].diff()**2 -
-                       2*l['l']*cos(l['dc'])*l['dc'].diff()*l['fd'].diff() -
-                       l['l']*sin(l['dc'])*l['fd'].diff().diff())*l['frame_d'].y +
-                      (l['r']*l['fd'].diff().diff() +
-                       l['l']*cos(l['dc'])*l['dc'].diff().diff() -
-                       l['l']*sin(l['dc'])*l['dc'].diff()**2 -
-                       l['l']*sin(l['dc'])*l['fd'].diff()**2)*l['frame_d'].z)
-    assert (a_e_f - expected_a_e_f).simplify() == 0
diff --git a/sympy/physics/continuum_mechanics/tests/test_beam.py b/sympy/physics/continuum_mechanics/tests/test_beam.py
index 17db45e4a7..48b8d2a773 100644
--- a/sympy/physics/continuum_mechanics/tests/test_beam.py
+++ b/sympy/physics/continuum_mechanics/tests/test_beam.py
@@ -472,10 +472,11 @@ def test_Beam3D():
 
     assert b.shear_force() == [0, -q*x, 0]
     assert b.bending_moment() == [0, 0, -m*x + q*x**2/2]
-    assert b.deflection() == [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l**2*q - 2*A*G*l*m
-            + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I)) + l*m*x**2/(4*E*I) - l*x**3*(A*G*l**2*q
-            - 2*A*G*l*m + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I)) - m*x**3/(6*E*I) + q*x**4/(24*E*I)
-            + l*x*(A*G*l**2*q - 2*A*G*l*m + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I))
+    assert b.deflection() == [0, -l**2*q*x**2/(12*E*I) + l**2*x**2*(A*G*l*(l*q - 2*m)
+            + 12*E*I*q)/(8*E*I*(A*G*l**2 + 12*E*I)) + l*m*x**2/(4*E*I)
+            - l*x**3*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(12*E*I*(A*G*l**2 + 12*E*I))
+            - m*x**3/(6*E*I) + q*x**4/(24*E*I)
+            + l*x*(A*G*l*(l*q - 2*m) + 12*E*I*q)/(2*A*G*(A*G*l**2 + 12*E*I))
             - q*x**2/(2*A*G), 0]
 
 
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 3785645add..bf45d2acc2 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1019,26 +1019,42 @@ def test_abs_invert_solvify():
 
 def test_linear_eq_to_matrix():
     x, y, z = symbols('x, y, z')
+    a, b, c, d, e, f, g, h, i, j, k, l = symbols('a:l')
+
     eqns1 = [2*x + y - 2*z - 3, x - y - z, x + y + 3*z - 12]
     eqns2 = [Eq(3*x + 2*y - z, 1), Eq(2*x - 2*y + 4*z, -2), -2*x + y - 2*z]
 
-    A, b = linear_eq_to_matrix(eqns1, x, y, z)
+    A, B = linear_eq_to_matrix(eqns1, x, y, z)
     assert A == Matrix([[2, 1, -2], [1, -1, -1], [1, 1, 3]])
-    assert b == Matrix([[3], [0], [12]])
+    assert B == Matrix([[3], [0], [12]])
 
-    A, b = linear_eq_to_matrix(eqns2, x, y, z)
+    A, B = linear_eq_to_matrix(eqns2, x, y, z)
     assert A == Matrix([[3, 2, -1], [2, -2, 4], [-2, 1, -2]])
-    assert b == Matrix([[1], [-2], [0]])
+    assert B == Matrix([[1], [-2], [0]])
 
     # Pure symbolic coefficients
-    from sympy.abc import a, b, c, d, e, f, g, h, i, j, k, l
-    eqns3 = [a*x + b*y + c*z - d, e*x + f*y + g*z - h, i*x + j*y + k*z - l]
+    eqns3 = [a*b*x + b*y + c*z - d, e*x + d*x + f*y + g*z - h, i*x + j*y + k*z - l]
     A, B = linear_eq_to_matrix(eqns3, x, y, z)
-    assert A == Matrix([[a, b, c], [e, f, g], [i, j, k]])
+    assert A == Matrix([[a*b, b, c], [d + e, f, g], [i, j, k]])
     assert B == Matrix([[d], [h], [l]])
 
-    # raise ValueError if no symbols are given
+    # raise ValueError if
+    # 1) no symbols are given
     raises(ValueError, lambda: linear_eq_to_matrix(eqns3))
+    # 2) there are duplicates
+    raises(ValueError, lambda: linear_eq_to_matrix(eqns3, [x, x, y]))
+    # 3) there are non-symbols
+    raises(ValueError, lambda: linear_eq_to_matrix(eqns3, [x, 1/a, y]))
+    # 4) a nonlinear term is detected in the original expression
+    raises(ValueError, lambda: linear_eq_to_matrix(Eq(1/x + x, 1/x)))
+
+    assert linear_eq_to_matrix(1, x) == (Matrix([[0]]), Matrix([[-1]]))
+    # issue 15195
+    assert linear_eq_to_matrix(x + y*(z*(3*x + 2) + 3), x) == (
+        Matrix([[3*y*z + 1]]), Matrix([[-y*(2*z + 3)]]))
+    assert linear_eq_to_matrix(Matrix(
+        [[a*x + b*y - 7], [5*x + 6*y - c]]), x, y) == (
+        Matrix([[a, b], [5, 6]]), Matrix([[7], [c]]))
 
 
 def test_linsolve():
@@ -1133,6 +1149,13 @@ def test_linsolve():
     Eqns = [8*kilo*newton + x + y, 28*kilo*newton*meter + 3*x*meter]
     assert linsolve(Eqns, x, y) == {(-28000*newton/3, 4000*newton/3)}
 
+    # linsolve fully expands expressions, so removable singularities
+    # and other nonlinearity does not raise an error
+    assert linsolve([Eq(x, x + y)], [x, y]) == {(x, 0)}
+    assert linsolve([Eq(1/x, 1/x + y)], [x, y]) == {(x, 0)}
+    assert linsolve([Eq(y/x, y/x + y)], [x, y]) == {(x, 0)}
+    assert linsolve([Eq(x*(x + 1), x**2 + y)], [x, y]) == {(y, y)}
+
 
 def test_solve_decomposition():
     x = Symbol('x')
@@ -1591,6 +1614,27 @@ def test_issue_10397():
     assert solveset(sqrt(x), x, S.Complexes) == FiniteSet(0)
 
 
+def test_issue_14987():
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [x**2], x))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [x*(-3/x + 1) + 2*y - a], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [(x**2 - 3*x)/(x - 3) - 3], x))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [(x + 1)**3 - x**3 - 3*x**2 + 7], x))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [x*(1/x + 1) + y], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [(x + 1)*y], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [Eq(1/x, 1/x + y)], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [Eq(y/x, y/x + y)], [x, y]))
+    raises(ValueError, lambda: linear_eq_to_matrix(
+        [Eq(x*(x + 1), x**2 + y)], [x, y]))
+
+
 def test_simplification():
     eq = x + (a - b)/(-2*a + 2*b)
     assert solveset(eq, x) == FiniteSet(S.Half)