Repository: sympy/sympy

Issue:
Assertion error using integrate in sympy 1.10.1, Python 3.10.5
I am getting this error

```
>python
Python 3.10.5 (main, Jun  6 2022, 18:49:26) [GCC 12.1.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from sympy import *
>>>F, G, H, a, b, c, d, e, f, g, h, i, k, m, n, p, r, s, t, x=symbols('F G H a b c d e f g h i k m n p r s t x')
>>>integrand=((2*x*exp(x)**2+(2*x**2-2)*exp(x)-2*x)*exp(exp(x))**2*ln(2*exp(x)+2*x)**2+((2*exp(x)*x+2*x)*exp(exp(x))**2+(-6*exp(x)**2*x**2+(-6*x**3+6*x)*exp(x)+6*x**2)*exp(exp(x)))*ln(2*exp(x)+2*x)+(-6*exp(x)*x**2-6*x**2)*exp(exp(x))+exp(x)*x**3+x**4)/(exp(x)*x**3+x**4)
>>>integrate(integrand,x)
```
gives

```
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 1566, in integrate
    return integral.doit(**doit_flags)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 613, in doit
    antideriv = self._eval_integral(
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 956, in _eval_integral
    result, i = risch_integrate(f, x, separate_integral=True,
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1831, in risch_integrate
    ans, i, b = integrate_primitive(fa, fd, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1460, in integrate_primitive
    q, i, b = integrate_primitive_polynomial(p, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1409, in integrate_primitive_polynomial
    rv = limited_integrate(aa, ad, [(Dta, Dtb)], DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 837, in limited_integrate
    h, A = param_rischDE(Fa, Fd, G, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 734, in param_rischDE
    h, B = param_poly_rischDE(a, b, r, n, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 559, in param_poly_rischDE
    return prde_cancel_liouvillian(b, q, n, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 481, in prde_cancel_liouvillian
    fi, Ai = param_rischDE(ba, bd, Qy, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 734, in param_rischDE
    h, B = param_poly_rischDE(a, b, r, n, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 559, in param_poly_rischDE
    return prde_cancel_liouvillian(b, q, n, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 481, in prde_cancel_liouvillian
    fi, Ai = param_rischDE(ba, bd, Qy, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 734, in param_rischDE
    h, B = param_poly_rischDE(a, b, r, n, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 544, in param_poly_rischDE
    return prde_no_cancel_b_large(b, q, n, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 355, in prde_no_cancel_b_large
    A = A.row_join(zeros(A.rows, m, DE.t)).col_join(c.row_join(-c))
  File "/usr/lib/python3.10/site-packages/sympy/polys/polymatrix.py", line 253, in col_join
    dm = DomainMatrix.vstack(self._dm, other._dm)
  File "/usr/lib/python3.10/site-packages/sympy/polys/matrices/domainmatrix.py", line 788, in vstack
    return DomainMatrix.from_rep(A.rep.vstack(*(Bk.rep for Bk in B)))
  File "/usr/lib/python3.10/site-packages/sympy/polys/matrices/sdm.py", line 794, in vstack
    assert Bkcols == cols
AssertionError

```

Versions

```
>which python
/usr/bin/python
>python --version
Python 3.10.5
>python
Python 3.10.5 (main, Jun  6 2022, 18:49:26) [GCC 12.1.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import sympy 
>>> sympy.__version__
'1.10.1'
>>> 


```

On Linux Manjaro


Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* integrals
  * A bug in the `risch` integration routine was fixed. Previously in certain conditions involving complicated integrands the routine would fail with `AssertionError`. ([#23959](https://github.com/sympy/sympy/pull/23959) by [@oscarbenjamin](https://github.com/oscarbenjamin))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.12.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->

    Fixes #23948

    #### Brief description of what is fixed or changed

    Ensure that while the matrix is still empty it has the right number of columns.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below between the BEGIN and END
    statements. The basic format is a bulleted list with the name of the subpackage
    and the release note for this PR. For example:

    * solvers
      * Added a new solver for logarithmic equations.

    * functions
      * Fixed a bug with log of integers.

    or if no release note(s) should be included use:

    NO ENTRY

    See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
    information on how to write release notes. The bot will check your release
    notes automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    * integrals
       * A bug in the `risch` integration routine was fixed. Previously in certain conditions involving complicated integrands the routine would fail with `AssertionError`.
    <!-- END RELEASE NOTES -->


</details><p>

asmeurer: This line is also probably wrong https://github.com/sympy/sympy/blob/88664e6e0b781d0a8b5347896af74b555e92891e/sympy/integrals/prde.py#L388
asmeurer: Also I'm surprised that `Matrix()` with no arguments works (it creates a 0x0 matrix). I guess that's been around for some time, but IMO it's not great. 
oscarbenjamin: > This line is also probably wrong

Can you think of a test-case that would hit that line?

> Also I'm surprised that `Matrix()` with no arguments works

Nothing like this surprises me any more.
github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [26f7bdbe]       [5eb59bda]
     <sympy-1.11^0>                 
-         965±5μs        614±0.9μs     0.64  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-        2.79±0ms         1.15±0ms     0.41  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.56±0.01ms         1.70±0ms     0.31  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

Patch:
diff --git a/doc/src/explanation/index.rst b/doc/src/explanation/index.rst
index 1306d26627..fb71649896 100644
--- a/doc/src/explanation/index.rst
+++ b/doc/src/explanation/index.rst
@@ -12,10 +12,9 @@ technical implementation details, and opinionated recommendations.
 **Content**
 
 .. toctree::
-   :maxdepth: 1
+   :maxdepth: 2
 
    gotchas.rst
-   solve_output.rst
    special_topics/index.rst
    active-deprecations.md
    glossary.md
diff --git a/doc/src/explanation/solve_output.rst b/doc/src/explanation/solve_output.rst
deleted file mode 100644
index 2815906739..0000000000
--- a/doc/src/explanation/solve_output.rst
+++ /dev/null
@@ -1,142 +0,0 @@
-.. _solve_output:
-
-====================
-Solve Output by Type
-====================
-
-The output of the :func:`~.solve` function can seem very unwieldy since it may appear to
-arbitrarily return one of six different types of output (in addition to raising
-errors). The reasons for this are historical and are biased toward human
-interaction rather than programmatic use. The type of output will depend on the
-type of equation(s) (and how they are entered) and the number of symbols that
-are provided (and how they are provided).
-
-    >>> from sympy import sqrt, exp, solve, Symbol, Eq
-    >>> from sympy.abc import x, y, z, a, b
-
-    The :func:`~.solve` function attempts to find all values for as many symbols as
-    possible that will make each expression given equal to zero. The format
-    of the output can be controlled by using the ``dict`` or ``set`` keyword:
-
-    >>> solve(x - 1, dict=True)
-    [{x: 1}]
-    >>> solve([x**2 - y, x + y - 6], set=True)
-    ([x, y], {(-3, 9), (2, 4)})
-
-    The following discussion provides an explanation for the output
-    obtained when not using those keywords.
-
-Empty List
-----------
-
-    When there is no solution, an empty list is returned.
-
-    >>> solve(sqrt(x) + 1)  # or solve(sqrt(x) + 1, dict=True)
-    []
-    >>> solve(sqrt(x) + 1, set=True)
-    ([x], set())
-
-List Of Values
---------------
-
-    A list of values is given when the symbol to solve for was
-    unambiguous in context because a) the equation was univariate or
-    b) a single symbol was specified as being of interest.
-
-    >>> solve(x**2 - 4)
-    [-2, 2]
-    >>> solve(x - y - 1, x)
-    [y + 1]
-
-Single Dictionary
------------------
-
-    A single dictionary with keys being symbols and values being the solutions
-    for those symbols is the result when equations are passed as a list and are
-    all linear in the symbols given. Note: such a system is automatically generated
-    for a single equation (not passed as a list) if there is an
-    undetermined-coefficients solution for the symbols specified. If this is not
-    what was intended, then pass the expression in a list.
-
-    >>> solve([x + y - 2, x - y + 2], x, y)
-    {x: 0, y: 2}
-    >>> eq = a*x - 2*x + b - 5
-    >>> solve(eq, {a, b})  # undetermined coefficients
-    {a: 2, b: 5}
-    >>> solve([eq], {a, b})  # algebraic
-    {a: -b/x + (2*x + 5)/x}
-
-List of Tuples
---------------
-
-    Each tuple in the list gives a solution for the symbols in the order
-    they were given. This format is used when a) a list of equations contains
-    at least one nonlinear equation or b) a list of symbols is given in a
-    well defined order. (This is also the format for the tuples in the set
-    returned when using the flag ``set=True``.)
-
-    >>> solve(x - 1, x, y)  # more than one symbol
-    [(1, y)]
-    >>> solve([x**2], x)  # list with nonlinear equation
-    [(0,)]
-    >>> solve([x**2 - 1], x)
-    [(-1,), (1,)]
-    >>> solve([x**2 - y, x - 3], x, y)  # nonlinear and multiple symbols
-    [(3, 9)]
-
-List of Dictionaries
---------------------
-
-    The list of dictionaries is returned when the expression was not
-    univariate or there was a nonlinear expression in a list *and* the order of
-    symbols would otherwise be ambiguous because a) no symbols were passed or
-    b) the symbols were passed as a set. (This is also the format selected with
-    ``dict=True``.)
-
-    >>> solve(x - y)
-    [{x: y}]
-    >>> solve([exp(x) - 1, x*(x - 1)])
-    [{x: 0}]
-    >>> system = [x + y - z, x**2 - y + z, exp(z) + 1/x + 1/y - 2]
-    >>> sol = solve(system[:2]); sol
-    [{x: -1, y: z + 1}, {x: 0, y: z}]
-
-    The dictionaries only contain values that are distinct from the keys.
-    In the last example above, there is no key for ``z`` in the dictionary
-    since only *two* of the three equations were insufficient to determine its
-    value. These solutions can be used to eliminate those variables from the
-    third equation, however, to give a relationship in a single variable that
-    can be solved (perhaps numerically) to obtain a full solution with the
-    advantage of only needing to guess a single value instead of three.
-
-        >>> from sympy import nsolve
-        >>> [system[-1].subs(s) for s in sol]
-        [exp(z) - 3 + 1/(z + 1), exp(z) + zoo + 1/z]
-        >>> z_eq = _[0]
-        >>> zsol = nsolve(z_eq, 1); zsol
-        0.906425478894557
-        >>> sol0 = {k: v.subs(z, zsol) for k, v in sol[0].items()}
-        >>> sol0[z] = zsol; sol0
-        {x: -1, y: 1.90642547889456, z: 0.906425478894557}
-
-Boolean or Relational
----------------------
-
-    A boolean expression is returned when a relational expression other
-    than an :class:`~.Equality` is given as an expression to solve. A single `Equality`
-    or a more complicated relational expression might be returned. The
-    use of :func:`~.solve` here is equivalent to passing the equation set and
-    symbols to :func:`~.reduce_inequalities` (and ``dict``, ``set``, and ``check``
-    flags are ignored).
-
-    >>> solve([x**2 > 4, x > 0])
-    (2 < x) & (x < oo)
-
-    >>> from sympy import Unequality as Ne
-    >>> solve([x**2 - 4, Ne(x, -2)])
-    Eq(x, 2)
-
-    Any returned `Equality` can be converted to a dictionary:
-
-    >>> {_.lhs: _.rhs}
-    {x: 2}
diff --git a/doc/src/guides/solving/index.md b/doc/src/guides/solving/index.md
index c50b9b9334..5a08094fb4 100644
--- a/doc/src/guides/solving/index.md
+++ b/doc/src/guides/solving/index.md
@@ -5,9 +5,6 @@ The Python package SymPy can symbolically solve equations, differential equation
 linear equations, nonlinear equations, matrix problems, inequalities,
 Diophantine equations, and evaluate integrals. SymPy can also solve numerically.
 
-The [solving guide](solving-guidance.md) provides suggestions for many
-types of solving tasks.
-
 Learn how to use SymPy computer algebra system to:
 
 | Description                                                  | Example                                                                                                                     | Solution |
@@ -34,6 +31,5 @@ in SymPy's terminology it would be ["evaluate an integral."](../../modules/integ
 ```{toctree}
 :hidden: true
 
-solving-guidance.md
 solve-equation-algebraically.md
 ```
diff --git a/doc/src/guides/solving/solve-equation-algebraically.md b/doc/src/guides/solving/solve-equation-algebraically.md
index 99a1fd792f..f0a7e6a4de 100644
--- a/doc/src/guides/solving/solve-equation-algebraically.md
+++ b/doc/src/guides/solving/solve-equation-algebraically.md
@@ -1,14 +1,15 @@
-# Solve an Equation Algebraically
+# Solve an equation algebraically
 
 Use SymPy to solve an equation algebraically (symbolically). For example,
 solving $x^2 = y$ for $x$ yields $x \in \{-\sqrt{y},\sqrt{y}\}$.
 
 Alternatives to consider:
-- SymPy can also [solve many other types of problems including sets of
-equations](index.md).
-- Some equations cannot be solved algebraically (either at all or by SymPy), so
-you may have to {func}`solve your equation numerically
-<sympy.solvers.solvers.nsolve>` instead.
+- SymPy can also
+[solve many other types of problems including sets of equations](index.md).
+- Some equations cannot be solved algebraically (either at all or by SymPy),
+so you may have to
+{func}`solve your equation numerically <sympy.solvers.solvers.nsolve>`
+instead.
 
 There are two high-level functions to solve equations, {func}`~.solve` and
 {func}`~.solveset`. Here is a simple example of each:
@@ -34,15 +35,15 @@ There are two high-level functions to solve equations, {func}`~.solve` and
 Here are recommendations on when to use:
 
 - {func}`~.solve`
-    - You want to get explicit symbolic representations of the different values
-    a variable could take that would satisfy the equation.
-    - You want to substitute those explicit solution values into other equations
-    or expressions involving the same variable using
+    - You want to get explicit symbolic representations of the different
+    values a variable could take that would satisfy the equation.
+    - You want to substitute those explicit solution values into other
+    equations or expressions involving the same variable using
     {meth}`~sympy.core.basic.Basic.subs`
 
 - {func}`~.solveset`
-    - You want to represent the solutions in a mathematically precise way, using
-    [mathematical sets](../../modules/sets.rst).
+    - You want to represent the solutions in a mathematically precise way,
+    using [mathematical sets](../../modules/sets.rst).
     - You want a representation of all the solutions, including if there are
     infinitely many.
     - You want a consistent input interface.
@@ -52,23 +53,53 @@ Here are recommendations on when to use:
 
 ## Guidance
 
-Refer to
-[](solving-guidance.md#include-the-variable-to-be-solved-for-in-the-function-call)
-and [](solving-guidance.md#ensure-consistent-formatting-from).
+### Include the variable to be solved for in the function call
+
+We recommend you include the variable to be solved for as the second argument
+for either function. While this is optional for equations with a single
+symbol, it is a good practice because it ensures
+SymPy will solve for the desired symbol. For example, you may expect the
+following to solve for $x$,
+and SymPy will solve for $y$:
+
+```py
+>>> from sympy.abc import x, y
+>>> from sympy import solve
+>>> solve(x**2 - y, dict=True)
+[{y: x**2}]
+```
+
+Specifying the variable to solve for ensures that SymPy solves for it:
+
+```py
+>>> from sympy.abc import x, y
+>>> from sympy import solve
+>>> solve(x**2 - y, x, dict=True)
+[{x: -sqrt(y)}, {x: sqrt(y)}]
+```
+
+### Ensure consistent formatting from {func}`~.solve` by using `dict=True`
 
-## Solve an Equation Algebraically
+{func}`~.solve` produces various output formats depending on the answer,
+unless you use `dict=True` to ensure the result will be formatted as a
+dictionary. We recommend using `dict=True`, especially if you want to
+extract information from the result programmatically.
 
-You can solve an equation in several ways. The examples below demonstrate using
-both {func}`~.solve` and {func}`~.solveset` where applicable. You can choose the
-function best suited to your equation.
+## Solve an equation using {func}`~.solve` or {func}`~.solveset`
 
-### Make Your Equation Into an Expression That Equals Zero
+You can solve an equation in several ways.
+The examples below demonstrate using both {func}`~.solve` and
+{func}`~.solveset` where applicable.
+You can choose the function best suited to your equation.
+
+### Make your equation into an expression that equals zero
 
 Use the fact that any expression not in an `Eq` (equation) is automatically
 assumed to equal zero (0) by the solving functions. You can rearrange the
-equation $x^2 = y$ to $x^2 - y = 0$, and solve that expression. This approach is
-convenient if you are interactively solving an expression which already equals
-zero, or an equation that you do not mind rearranging to $expression = 0$.
+equation $x^2 = y$ to $x^2 - y = 0$, and solve that expression. This approach
+is convenient if you are interactively solving an expression which already
+equals zero, or an equation that you do not mind rearranging to
+$expression = 0$.
 
 ```py
 >>> from sympy import solve, solveset
@@ -79,12 +110,12 @@ zero, or an equation that you do not mind rearranging to $expression = 0$.
 {-sqrt(y), sqrt(y)}
 ```
 
-### Put Your Equation Into `Eq` Form
+### Put your equation into `Eq` form
 
-Put your equation into `Eq` form, then solve the `Eq`. This approach is
-convenient if you are interactively solving an equation which you already have
-in the form of an equation, or which you think of as an equality. It also helps
-to prevent sign errors when subtracting one side from the other.
+Put your equation into `Eq` form, then solve the `Eq`.
+This approach is convenient if you are interactively solving an equation
+which you already have in the form of an equation,
+or which you think of as an equality.
 
 ```py
 >>> from sympy import Eq, solve, solveset
@@ -104,7 +135,7 @@ sqrt(y)
 -sqrt(y)
 ```
 
-### Restrict the Domain of Solutions
+### Restrict the domain of solutions
 
 By default, SymPy will return solutions in the complex domain, which also
 includes purely real and imaginary values. Here, the first two solutions are
@@ -119,8 +150,8 @@ real, and the last two are imaginary:
 {-4, 4, -4*I, 4*I}
 ```
 
-To restrict returned solutions to real numbers, or another domain or range, the
-different solving functions use different methods.
+To restrict returned solutions to real numbers, or another domain or range,
+the different solving functions use different methods.
 
 For {func}`~.solve`, place an assumption on the symbol to be solved for, $x$
 
@@ -146,8 +177,8 @@ such as a list comprehension:
 [1, 4]
 ```
 
-For {func}`~.solveset`, limit the output domain in the function call by setting
-a domain
+For {func}`~.solveset`, limit the output domain in the function call by
+setting a domain
 
 ```py
 >>> from sympy import S, solveset
@@ -177,7 +208,7 @@ and if you restrict the solutions to a domain in which there are no solutions,
 EmptySet
 ```
 
-### Explicitly Represent Infinite Sets of Possible Solutions
+### Explicitly represent infinite sets of possible solutions using {func}`~.solveset`
 
 {func}`~.solveset` [can represent infinite sets of possible
 solutions](why-solveset) and express them in standard mathematical notation, for
@@ -204,24 +235,25 @@ However, {func}`~.solve` will return only a finite number of solutions:
 2*pi
 ```
 
-{func}`~.solve` tries to return just enough solutions so that all (infinitely
-many) solutions can generated from the returned solutions by adding integer
-multiples of the {func}`~.periodicity` of the equation, here $2\pi$.
+{func}`~.solve` tries to return just enough solutions so that all
+(infinitely many) solutions can generated from the returned solutions
+by adding integer multiples of the {func}`~.periodicity` of the
+equation, here $2\pi$.
 
-## Use the Solution Result
+## Use the solution result
 
-### Substitute Solutions From {func}`~.solve` Into an Expression
+### Substitute solutions from {func}`~.solve` into an expression
 
 You can substitute solutions from {func}`~.solve` into an expression.
 
-A common use case is finding the critical points and values for a function $f$.
-At the critical points, the {class}`~.Derivative` equals zero (or is undefined).
-You can then obtain the function values at those critical points by substituting
-the critical points back into the function using
-{meth}`~sympy.core.basic.Basic.subs`. You can also tell if the critical point is
-a maxima or minima by substituting the values into the expression for the second
-derivative: a negative value indicates a maximum, and a positive value indicates
-a minimum.
+A common use case is finding the critical points and values for a function
+$f$. At the critical points, the {class}`~.Derivative` equals zero (or is
+undefined). You can then obtain the function values at those critical points
+by substituting the critical points back into the function using
+{meth}`~sympy.core.basic.Basic.subs`. You can also tell if the critical point
+is a maxima or minima by substituting the values into the expression for the
+second derivative: a negative value indicates a maximum, and a positive value
+indicates a minimum.
 
 ```py
 >>> from sympy.abc import x
@@ -243,7 +275,7 @@ a minimum.
 4
 ```
 
-### {func}`~.solveset` Solution Sets Cannot Necessarily Be Interrogated Programmatically
+### {func}`~.solveset` solution sets cannot necessarily be interrogated programmatically
 
 If {func}`~.solveset` returns a finite set (class {class}`~.FiniteSet`), you can
 iterate through the solutions:
@@ -259,7 +291,8 @@ iterate through the solutions:
 [sqrt(y), -sqrt(y)]
 ```
 
-However, for more complex results, it may not be possible to list the solutions:
+However, for more complex results, it may not be possible to list the
+solutions:
 
 ```py
 >>> from sympy import S, solveset, symbols
@@ -308,36 +341,70 @@ Intersection({-sqrt(y), sqrt(y)}, Reals)
 [sqrt(y), -sqrt(y)]
 ```
 
-## Options That Can Speed up {func}`~.solve`
+## Options that can speed up {func}`~.solve`
+
+### Include solutions making any denominator zero by using `check=False`
+
+Normally, {func}`~.solve` checks whether any solutions make any denominator
+zero, and automatically excludes them. If you want to include those
+solutions, and speed up {func}`~.solve` (at the risk of obtaining invalid
+solutions), set `check=False`:
+
+```py
+>>> from sympy import Symbol, sin, solve
+>>> x = Symbol("x")
+>>> solve(sin(x)/x)  # 0 is excluded
+[pi]
+>>> solve(sin(x)/x, check=False)  # 0 is not excluded
+[0, pi]
+```
+
+### Do not simplify solutions by using `simplify=False`
 
-Refer to [solving guidance](solving-guidance.md#options-that-can-speed-up).
+Normally, {func}`~.solve` simplifies all but polynomials of order 3 or
+greater before returning them and (if `check` is not False) uses the general
+{func}`simplify <sympy.simplify.simplify.simplify>` function on the solutions
+and the expression obtained when they are substituted into the function which
+should be zero. If you do not want the solutions simplified, and want to
+speed up {func}`~.solve`, use `simplify=False`.
 
-## Not All Equations Can Be Solved
+```py
+>>> from sympy import solve
+>>> from sympy.abc import x, y
+>>> expr = x**2 - (y**5 - 3*y**3 + y**2 - 3)
+>>> solve(expr, x, dict=True)
+[{x: -sqrt(y**5 - 3*y**3 + y**2 - 3)}, {x: sqrt(y**5 - 3*y**3 + y**2 - 3)}]
+>>> solve(expr, x, dict=True, simplify=False)
+[{x: -sqrt((y + 1)*(y**2 - 3)*(y**2 - y + 1))}, {x: sqrt((y + 1)*(y**2 - 3)*(y**2 - y + 1))}]
+```
 
-### Equations With No Closed-Form Solution
+## Not all equations can be solved
 
-Some equations have no closed-form solution, in which case SymPy may return an
-empty set or give an error. For example, the following [transcendental
-equation](https://en.wikipedia.org/wiki/Transcendental_equation) has no
-closed-form solution:
+### Equations with no solution
+
+Some equations have no solution, in which case SymPy may return an empty set.
+For example, the equation $x - 7 = x + 2$ reduces to $-7 = 2$, which has no
+solution because no value of $x$ will make it true:
 
 ```py
->>> from sympy import cos, solve
+>>> from sympy import solve, Eq
 >>> from sympy.abc import x
->>> solve(cos(x) - x, x, dict=True)
-Traceback (most recent call last):
-    ...
-NotImplementedError: multiple generators [x, cos(x)]
-No algorithms are implemented to solve equation -x + cos(x)
+>>> eqn = Eq(x - 7, x + 2)
+>>> solve(eqn, x)
+[]
 ```
 
-### Equations Which Have a Closed-Form Solution, and SymPy Cannot Solve
+So if SymPy returns an empty list, you may want to check whether there is a
+mistake in the equation.
+
+### Equations which have an analytical solution, and SymPy cannot solve
 
-It is also possible that there is an algebraic solution to your equation, and
-SymPy has not implemented an appropriate algorithm. If that happens, or SymPy
-returns an empty set or list when there is a mathematical solution (indicating a
-bug in SymPy), please post it on the [mailing
-list](https://groups.google.com/g/sympy), or open an issue on [SymPy's GitHub
-page](https://github.com/sympy/sympy/issues). Until the issue is resolved, you
-can {func}`solve your equation numerically <sympy.solvers.solvers.nsolve>`
+It is also possible that there is an algebraic solution to your equation,
+and SymPy has not implemented an appropriate algorithm.
+If that happens, or SymPy returns an empty set or list when there is a
+mathematical solution (indicating a bug in SymPy), please post it on the
+[mailing list](https://groups.google.com/g/sympy), or open an issue on
+[SymPy's GitHub page](https://github.com/sympy/sympy/issues). Until the issue
+is resolved, you can
+{func}`solve your equation numerically <sympy.solvers.solvers.nsolve>`
 instead.
diff --git a/doc/src/guides/solving/solving-guidance.md b/doc/src/guides/solving/solving-guidance.md
deleted file mode 100644
index 9991a9cbe4..0000000000
--- a/doc/src/guides/solving/solving-guidance.md
+++ /dev/null
@@ -1,261 +0,0 @@
-# Solving Guidance
-
-These guidelines apply to many types of solving.
-
-## Numeric Solutions
-
-### Equations With no Closed-Form Solution
-
-The vast majority of arbitrary nonlinear equations have no closed-form solution.
-The classes of equations that are solvable are basically:
-1. Linear equations
-2. Polynomials, except where limited by the [Abel-Ruffini
-   theorem](https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem) (learn
-   more about solving polynomials using a {class}`~.GroebnerBasis`)
-3. Equations that can be solved by inverting some transcendental functions
-4. Problems that can be transformed into the cases above (e.g., by turning
-trigonometric functions into polynomials)
-5. A few other special cases that can be solved with something like the
-{class}`Lambert W function <sympy.functions.elementary.exponential.LambertW>`
-6. Equations that you can {func}`~sympy.polys.polytools.decompose` via any of
-   the above
-
-SymPy may reflect that your equation has no solutions that can be expressed
-algebraically (symbolically), or that SymPy lacks an algorithm to find a
-closed-form solution that does exist, by returning an error such as
-`NotImplementedError`:
-
-```py
->>> from sympy import solve, cos
->>> from sympy.abc import x
->>> solve(cos(x) - x, x, dict=True)
-Traceback (most recent call last):
-  ...
-NotImplementedError: multiple generators [x, cos(x)]
-No algorithms are implemented to solve equation -x + cos(x)
-```
-
-so you may have to solve your equation numerically instead, for example using
-{func}`~.nsolve`
-
-```py
->>> from sympy import nsolve, cos
->>> from sympy.abc import x
->>> nsolve(cos(x) - x, x, 2)
-0.739085133215161
-```
-
-If you receive non-closed-form solutions such as {class}`CRootOf()
-<sympy.polys.rootoftools.ComplexRootOf>` (which represents an indexed complex
-root of a polynomial), you can evaluate them numerically using
-{func}`~sympy.core.evalf`:
-
-```py
->>> from sympy import solve
->>> from sympy.abc import x
->>> solutions = solve(x**5 - x - 1, x, dict=True)
->>> solutions
-[{x: CRootOf(x**5 - x - 1, 0)}, {x: CRootOf(x**5 - x - 1, 1)}, {x: CRootOf(x**5 - x - 1, 2)}, {x: CRootOf(x**5 - x - 1, 3)}, {x: CRootOf(x**5 - x - 1, 4)}]
->>> [solution[x].evalf(3) for solution in solutions]
-[1.17, -0.765 - 0.352*I, -0.765 + 0.352*I, 0.181 - 1.08*I, 0.181 + 1.08*I]
-```
-
-### When You Might Prefer a Numeric Solution
-Even if your problem has a closed-form solution, you might prefer a numeric
-solution.
-
-Solving functions such as {func}`~.solve` and {func}`~.solveset` will not try to
-find a numeric solution, only a mathematically-exact symbolic solution. So if
-you want a numeric solution, consider {func}`~.nsolve`.
-
-In some situations, even though a closed-form solution is available, it may be
-too cumbersome to be desirable. In that case, you can use
-{func}`~sympy.core.evalf` instead if a numerical solution is acceptable. For
-example, the following solution set contains more than 40 terms total when
-expressed exactly (scroll horizontally in the code block below if you want to
-view them all), compared to eight when expressed numerically:
-
-```py
->>> from sympy import symbols, solve
->>> x = symbols('x')
->>> solutions = solve(x**4 + 10*x**2 + x + 1, x, dict=True)
->>> solutions
-[{x: -sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2 - sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) + 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2}, {x: sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2 - sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) - 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2}, {x: sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) - 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2 + sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2}, {x: sqrt(-40/3 - 2*(1307/432 + sqrt(434607)*I/144)**(1/3) + 2/sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3)) - 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)))/2 - sqrt(-20/3 + 56/(9*(1307/432 + sqrt(434607)*I/144)**(1/3)) + 2*(1307/432 + sqrt(434607)*I/144)**(1/3))/2}]
->>> for solution in solutions:
-...     solution[x].evalf()
--0.0509758447494279 + 0.313552108895239*I
-0.0509758447494279 + 3.14751999969868*I
-0.0509758447494279 - 3.14751999969868*I
--0.0509758447494279 - 0.313552108895239*I
-```
-
-In other situations, even if the exact solution has few terms, you may want a
-numeric solution so you know its approximate numerical value. For example, it
-may be difficult to estimate that $\sqrt{2} e^{\pi}/2$ is approximately $16$:
-
-```py
->>> from sympy import pi, sqrt, exp, solve, evalf
->>> shorter = solve(sqrt(2)*x - exp(pi), x, dict=True)
->>> shorter
-[{x: sqrt(2)*exp(pi)/2}]
->>> [solution[x].evalf(3) for solution in shorter]
-[16.4]
-```
-
-## Use Exact Values
-
-If you want to preserve the exact mathematical values of symbols such as
-transcendental numbers and [square roots](symbolic-computation), define them so
-that SymPy can interpret them symbolically, for example use SymPy's
-{class}`.Pi`:
-
-```py
->>> from sympy import symbols, solve, pi
->>> x = symbols('x')
->>> solve(x**2 - pi, x, dict=True)
-[{x: -sqrt(pi)}, {x: sqrt(pi)}]
-```
-
-If you use the standard Python math version of $\pi$, Python will pass that
-inexact value to SymPy, leading to an inexact, numerical solution:
-
-```py
->>> from sympy import symbols, solve
->>> from math import pi
->>> x = symbols('x')
->>> solve(x**2 - pi, x, dict=True)
-[{x: -1.77245385090552}, {x: 1.77245385090552}]
-```
-
-To use exact values for numbers such as $6.2$ or $1/2$, refer to
-[](python-vs-sympy-numbers).
-
-In certain cases, using an inexact value will prevent SymPy from finding a
-result. For example, this exact equation can be solved:
-
-```py
->>> from sympy import symbols, solve, sqrt
->>> x = symbols('x')
->>> eq = x**sqrt(2) - 2
->>> solve(eq, x, dict=True)
-[{x: 2**(sqrt(2)/2)}]
-```
-
-but if you use the inexact equation `eq = x**1.4142135623730951 - 2`, SymPy will
-not return a result despite attempting for a long time. 
-
-## Include the Variable to be Solved for in the Function Call
-
-We recommend you include the variable to be solved for as the second argument
-for solving functions including {func}`~.solve` and {func}`~.solveset`. While
-this is optional for univariate equations, it is a good practice because it
-ensures SymPy will solve for the desired symbol. For example, you might be
-interested in a solution for $x$, but SymPy solves for $y$:
-
-```py
->>> from sympy.abc import x, y
->>> from sympy import solve
->>> solve(x**2 - y, dict=True)
-[{y: x**2}]
-```
-
-Specifying the variable to solve for ensures that SymPy solves for it:
-
-```py
->>> from sympy.abc import x, y
->>> from sympy import solve
->>> solve(x**2 - y, x, dict=True)
-[{x: -sqrt(y)}, {x: sqrt(y)}]
-```
-
-## Ensure Consistent Formatting From {func}`~.solve`
-
-{func}`~.solve` produces a variety of output as explained in
-{ref}`solve_output`. Using `dict=True` will give a consistent output format
-which is especially important when extracting information about the solution
-programmatically.
-
-To extract the solutions, you can iterate through the list of dictionaries:  
-    
-```py
->>> from sympy import parse_expr, solve, solveset
->>> from sympy.abc import x
->>> expr = "x^2 = y"
->>> parsed = parse_expr(expr, transformations="all")
->>> parsed
-Eq(x**2, y)
->>> solutions = solve(parsed, x, dict=True)
->>> [solution[x] for solution in solutions]
-[-sqrt(y), sqrt(y)]
->>> solveset(parsed, x)
-{-sqrt(y), sqrt(y)}
-```
-
-## Options That Can Speed up {func}`~.solve`
-
-### Include Solutions Making Any Denominator Zero
-
-Normally, {func}`~.solve` checks whether any solutions make any denominator
-zero, and automatically excludes them. If you want to include those solutions,
-and speed up {func}`~.solve` (at the risk of obtaining invalid solutions), set
-`check=False`:
-
-```py
->>> from sympy import Symbol, sin, solve
->>> x = Symbol("x")
->>> solve(sin(x)/x, x, dict=True) # 0 is excluded
-[{x: pi}]
->>> solve(sin(x)/x, x, dict=True, check=False) # 0 is not excluded
-[{x: 0}, {x: pi}]
-```
-
-### Do Not Simplify Solutions
-
-Normally, {func}`~.solve` simplifies many results before returning them and (if
-`check` is not False) uses the general {func}`~sympy.simplify.simplify.simplify`
-function on the solutions and the expression obtained when they are substituted
-into the function which should be zero. If you do not want the solutions
-simplified, and want to speed up {func}`~.solve`, use `simplify=False`.
-
-```py
->>> from sympy import solve
->>> from sympy.abc import x, y
->>> expr = x**2 - (y**5 - 3*y**3 + y**2 - 3)
->>> solve(expr, x, dict=True)
-[{x: -sqrt(y**5 - 3*y**3 + y**2 - 3)}, {x: sqrt(y**5 - 3*y**3 + y**2 - 3)}]
->>> solve(expr, x, dict=True, simplify=False)
-[{x: -sqrt((y + 1)*(y**2 - 3)*(y**2 - y + 1))}, {x: sqrt((y + 1)*(y**2 - 3)*(y**2 - y + 1))}]
-```
-
-## Parse a String Representing the Equation
-
-If you are creating the expression yourself, we advise [against using string
-parsing to create expressions](
-https://github.com/sympy/sympy/wiki/Idioms-and-Antipatterns#strings-as-input).
-But if you are programmatically reading in a string, this approach is
-convenient.
-
-You can parse a string representing the equation into a form that SymPy can
-understand (for example, {class}`~sympy.core.relational.Eq` form), then solve
-the parsed expression. Parsing an equation from a string requires you to use
-{func}`transformations <sympy.parsing.sympy_parser.parse_expr>` for SymPy to
-- interpret equals signs
-- create symbols from your variables
--  use more mathematical (rather than standard Python) notation, for example the
-exponent operator can be parsed from `^` rather than having to use Python's
-`**`.
-
-If you already have the equation in {class}`~sympy.core.relational.Eq`
-(equation) form, you can parse that string:
-
-```py
->>> from sympy import parse_expr, solve, solveset
->>> from sympy.abc import x
->>> expr = "Eq(x^2, y)"
->>> parsed = parse_expr(expr, transformations="all")
->>> parsed
-Eq(x**2, y)
-```
-
-SymPy can also parse [LaTeX](https://www.latex-project.org/) into expressions
-using {func}`~.parse_latex`. 
\ No newline at end of file
diff --git a/doc/src/tutorials/intro-tutorial/intro.rst b/doc/src/tutorials/intro-tutorial/intro.rst
index e095d83824..7471a0e3c0 100644
--- a/doc/src/tutorials/intro-tutorial/intro.rst
+++ b/doc/src/tutorials/intro-tutorial/intro.rst
@@ -2,8 +2,6 @@
  Introduction
 ==============
 
-.. _symbolic-computation:
-
 What is Symbolic Computation?
 =============================
 
diff --git a/sympy/core/sympify.py b/sympy/core/sympify.py
index b86897cc32..5364873550 100644
--- a/sympy/core/sympify.py
+++ b/sympy/core/sympify.py
@@ -190,7 +190,7 @@ def sympify(a, locals=None, convert_xor=True, strict=False, rational=False,
     multi-letter names that are defined in ``abc``).
 
     >>> from sympy.abc import _clash1
-    >>> set(_clash1)  # if this fails, see issue #23903
+    >>> set(_clash1)
     {'E', 'I', 'N', 'O', 'Q', 'S'}
     >>> sympify('I & Q', _clash1)
     I & Q
diff --git a/sympy/integrals/prde.py b/sympy/integrals/prde.py
index d8182155bf..80a8b44e23 100644
--- a/sympy/integrals/prde.py
+++ b/sympy/integrals/prde.py
@@ -344,10 +344,9 @@ def prde_no_cancel_b_large(b, Q, n, DE):
 
     if all(qi.is_zero for qi in Q):
         dc = -1
-        M = zeros(0, 2, DE.t)
     else:
         dc = max([qi.degree(DE.t) for qi in Q])
-        M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i), DE.t)
+    M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i), DE.t)
     A, u = constant_system(M, zeros(dc + 1, 1, DE.t), DE)
     c = eye(m, DE.t)
     A = A.row_join(zeros(A.rows, m, DE.t)).col_join(c.row_join(-c))
@@ -385,10 +384,9 @@ def prde_no_cancel_b_small(b, Q, n, DE):
             Q[i] = Q[i] - derivation(si, DE) - b*si
         if all(qi.is_zero for qi in Q):
             dc = -1
-            M = Matrix()
         else:
             dc = max([qi.degree(DE.t) for qi in Q])
-            M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i), DE.t)
+        M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i), DE.t)
         A, u = constant_system(M, zeros(dc + 1, 1, DE.t), DE)
         c = eye(m, DE.t)
         A = A.row_join(zeros(A.rows, m, DE.t)).col_join(c.row_join(-c))
diff --git a/sympy/physics/continuum_mechanics/beam.py b/sympy/physics/continuum_mechanics/beam.py
index 87dd9a6a80..445ee54257 100644
--- a/sympy/physics/continuum_mechanics/beam.py
+++ b/sympy/physics/continuum_mechanics/beam.py
@@ -2389,11 +2389,9 @@ def __init__(self, length, elastic_modulus, shear_modulus, second_moment,
         self.area = area
         self._load_vector = [0, 0, 0]
         self._moment_load_vector = [0, 0, 0]
-        self._torsion_moment = {}
         self._load_Singularity = [0, 0, 0]
         self._slope = [0, 0, 0]
         self._deflection = [0, 0, 0]
-        self._angular_deflection = 0
 
     @property
     def shear_modulus(self):
@@ -2558,11 +2556,6 @@ def apply_moment_load(self, value, start, order, dir="y"):
         if dir == "x":
             if not order == -2:
                 self._moment_load_vector[0] += value
-            else:
-                if start in list(self._torsion_moment):
-                    self._torsion_moment[start] += value
-                else:
-                    self._torsion_moment[start] = value
             self._load_Singularity[0] += value*SingularityFunction(x, start, order)
         elif dir == "y":
             if not order == -2:
@@ -2679,42 +2672,6 @@ def torsional_moment(self):
         """
         return self.bending_moment()[0]
 
-    def solve_for_torsion(self):
-        """
-        Solves for the angular deflection due to the torsional effects of
-        moments being applied in the x-direction i.e. out of or into the beam.
-
-        Here, a positive torque means the direction of the torque is positive
-        i.e. out of the beam along the beam-axis. Likewise, a negative torque
-        signifies a torque into the beam cross-section.
-
-        Examples
-        ========
-
-        >>> from sympy.physics.continuum_mechanics.beam import Beam3D
-        >>> from sympy import symbols
-        >>> l, E, G, I, A, x = symbols('l, E, G, I, A, x')
-        >>> b = Beam3D(20, E, G, I, A, x)
-        >>> b.apply_moment_load(4, 4, -2, dir='x')
-        >>> b.apply_moment_load(4, 8, -2, dir='x')
-        >>> b.apply_moment_load(4, 8, -2, dir='x')
-        >>> b.solve_for_torsion()
-        >>> b.angular_deflection().subs(x, 3)
-        18/(G*I)
-        """
-        x = self.variable
-        sum_moments = 0
-        for point in list(self._torsion_moment):
-            sum_moments += self._torsion_moment[point]
-        list(self._torsion_moment).sort()
-        pointsList = list(self._torsion_moment)
-        torque_diagram = Piecewise((sum_moments, x<=pointsList[0]), (0, x>=pointsList[0]))
-        for i in range(len(pointsList))[1:]:
-            sum_moments -= self._torsion_moment[pointsList[i-1]]
-            torque_diagram += Piecewise((0, x<=pointsList[i-1]), (sum_moments, x<=pointsList[i]), (0, x>=pointsList[i]))
-        integrated_torque_diagram = integrate(torque_diagram)
-        self._angular_deflection =  integrated_torque_diagram/(self.shear_modulus*self.polar_moment())
-
     def solve_slope_deflection(self):
         x = self.variable
         l = self.length
@@ -2797,13 +2754,6 @@ def deflection(self):
         """
         return self._deflection
 
-    def angular_deflection(self):
-        """
-        Returns a function in x depicting how the angular deflection, due to moments
-        in the x-axis on the beam, varies with x.
-        """
-        return self._angular_deflection
-
     def _plot_shear_force(self, dir, subs=None):
 
         shear_force = self.shear_force()
diff --git a/sympy/printing/aesaracode.py b/sympy/printing/aesaracode.py
index 6c4049b27a..8ac029d892 100644
--- a/sympy/printing/aesaracode.py
+++ b/sympy/printing/aesaracode.py
@@ -19,7 +19,7 @@
     mapping = {
             sympy.Add: aet.add,
             sympy.Mul: aet.mul,
-            sympy.Abs: aet.abs,
+            sympy.Abs: aet.abs_,
             sympy.sign: aet.sgn,
             sympy.ceiling: aet.ceil,
             sympy.floor: aet.floor,
@@ -51,10 +51,10 @@
             sympy.StrictLessThan: aet.lt,
             sympy.LessThan: aet.le,
             sympy.GreaterThan: aet.ge,
-            sympy.And: aet.bitwise_and,  # bitwise
-            sympy.Or: aet.bitwise_or,  # bitwise
+            sympy.And: aet.and_,  # bitwise
+            sympy.Or: aet.or_,  # bitwise
             sympy.Not: aet.invert,  # bitwise
-            sympy.Xor: aet.bitwise_xor,  # bitwise
+            sympy.Xor: aet.xor,  # bitwise
             sympy.Max: aet.maximum,  # Sympy accept >2 inputs, Aesara only 2
             sympy.Min: aet.minimum,  # Sympy accept >2 inputs, Aesara only 2
             sympy.conjugate: aet.conj,
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 321ce52604..99afd1be4d 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -382,107 +382,77 @@ def solve(f, *symbols, **flags):
         - piecewise combinations of the above
         - systems of linear and polynomial equations
         - systems containing relational expressions
-        - systems implied by undetermined coefficients
 
     Examples
     ========
 
-    The default output varies according to the input and might
-    be a list (possibly empty), a dictionary, a list of
-    dictionaries or tuples, or an expression involving relationals.
-    For specifics regarding different forms of output that may appear, see :ref:`solve_output`.
-    Let it suffice here to say that to obtain a uniform output from
-    `solve` use ``dict=True`` or ``set=True`` (see below).
+    The output varies according to the input and can be seen by example:
 
-        >>> from sympy import solve, Poly, Eq, Matrix, Symbol
+        >>> from sympy import solve, Poly, Eq, Function, exp
         >>> from sympy.abc import x, y, z, a, b
+        >>> f = Function('f')
+
+    Boolean or univariate Relational:
+
+        >>> solve(x < 3)
+        (-oo < x) & (x < 3)
+
+
+    To always get a list of solution mappings, use flag dict=True:
+
+        >>> solve(x - 3, dict=True)
+        [{x: 3}]
+        >>> sol = solve([x - 3, y - 1], dict=True)
+        >>> sol
+        [{x: 3, y: 1}]
+        >>> sol[0][x]
+        3
+        >>> sol[0][y]
+        1
+
 
-    The expressions that are passed can be Expr, Equality, or Poly
-    classes (or lists of the same); a Matrix is considered to be a
-    list of all the elements of the matrix:
+    To get a list of *symbols* and set of solution(s) use flag set=True:
 
+        >>> solve([x**2 - 3, y - 1], set=True)
+        ([x, y], {(-sqrt(3), 1), (sqrt(3), 1)})
+
+
+    Single expression and single symbol that is in the expression:
+
+        >>> solve(x - y, x)
+        [y]
         >>> solve(x - 3, x)
         [3]
         >>> solve(Eq(x, 3), x)
         [3]
         >>> solve(Poly(x - 3), x)
         [3]
-        >>> solve(Matrix([[x, x + y]]), x, y) == solve([x, x + y], x, y)
-        True
-
-    If no symbols are indicated to be of interest and the equation is
-    univariate, a list of values is returned; otherwise, the keys in
-    a dictionary will indicate which (of all the variables used in
-    the expression(s)) variables and solutions were found:
-
-        >>> solve(x**2 - 4)
-        [-2, 2]
-        >>> solve((x - a)*(y - b))
-        [{a: x}, {b: y}]
-        >>> solve([x - 3, y - 1])
-        {x: 3, y: 1}
-        >>> solve([x - 3, y**2 - 1])
-        [{x: 3, y: -1}, {x: 3, y: 1}]
-
-    If you pass symbols for which solutions are sought, the output will vary
-    depending on the number of symbols you passed, whether you are passing
-    a list of expressions or not, and whether a linear system was solved.
-    Uniform output is attained by using ``dict=True`` or ``set=True``.
-
-        >>> #### *** feel free to skip to the stars below *** ####
-        >>> from sympy import TableForm
-        >>> h = [None, ';|;'.join(['e', 's', 'solve(e, s)', 'solve(e, s, dict=True)',
-        ... 'solve(e, s, set=True)']).split(';')]
-        >>> t = []
-        >>> for e, s in [
-        ...         (x - y, y),
-        ...         (x - y, [x, y]),
-        ...         (x**2 - y, [x, y]),
-        ...         ([x - 3, y -1], [x, y]),
-        ...         ]:
-        ...     how = [{}, dict(dict=True), dict(set=True)]
-        ...     res = [solve(e, s, **f) for f in how]
-        ...     t.append([e, '|', s, '|'] + [res[0], '|', res[1], '|', res[2]])
-        ...
-        >>> # ******************************************************* #
-        >>> TableForm(t, headings=h, alignments="<")
-        e              | s      | solve(e, s)  | solve(e, s, dict=True) | solve(e, s, set=True)
-        ---------------------------------------------------------------------------------------
-        x - y          | y      | [x]          | [{y: x}]               | ([y], {(x,)})
-        x - y          | [x, y] | [(y, y)]     | [{x: y}]               | ([x, y], {(y, y)})
-        x**2 - y       | [x, y] | [(x, x**2)]  | [{y: x**2}]            | ([x, y], {(x, x**2)})
-        [x - 3, y - 1] | [x, y] | {x: 3, y: 1} | [{x: 3, y: 1}]         | ([x, y], {(3, 1)})
+        >>> solve(x**2 - y**2, x, set=True)
+        ([x], {(-y,), (y,)})
+        >>> solve(x**4 - 1, x, set=True)
+        ([x], {(-1,), (1,), (-I,), (I,)})
 
-        * If any equation does not depend on the symbol(s) given, it will be
-          eliminated from the equation set and an answer may be given
-          implicitly in terms of variables that were not of interest:
+    Single expression with no symbol that is in the expression:
 
-            >>> solve([x - y, y - 3], x)
-            {x: y}
-
-    When you pass all but one of the free symbols, an attempt
-    is made to find a single solution based on the method of
-    undetermined coefficients. If it succeeds, a dictionary of values
-    is returned. If you want an algebraic solutions for one
-    or more of the symbols, pass the expression to be solved in a list:
-
-        >>> e = a*x + b - 2*x - 3
-        >>> solve(e, [a, b])
-        {a: 2, b: 3}
-        >>> solve([e], [a, b])
-        {a: -b/x + (2*x + 3)/x}
-
-    When there is no solution for any given symbol which will make all
-    expressions zero, the empty list is returned (or an empty set in
-    the tuple when ``set=True``):
-
-        >>> from sympy import sqrt
         >>> solve(3, x)
         []
         >>> solve(x - 3, y)
         []
-        >>> solve(sqrt(x) + 1, x, set=True)
-        ([x], set())
+
+    Single expression with no symbol given. In this case, all free *symbols*
+    will be selected as potential *symbols* to solve for. If the equation is
+    univariate then a list of solutions is returned; otherwise - as is the case
+    when *symbols* are given as an iterable of length greater than 1 - a list of
+    mappings will be returned:
+
+        >>> solve(x - 3)
+        [3]
+        >>> solve(x**2 - y**2)
+        [{x: -y}, {x: y}]
+        >>> solve(z**2*x**2 - z**2*y**2)
+        [{x: -y}, {x: y}, {z: 0}]
+        >>> solve(z**2*x - z**2*y**2)
+        [{x: y**2}, {z: 0}]
 
     When an object other than a Symbol is given as a symbol, it is
     isolated algebraically and an implicit solution may be obtained.
@@ -491,9 +461,6 @@ def solve(f, *symbols, **flags):
     work if the specified object can be replaced with a Symbol using the
     subs method:
 
-        >>> from sympy import exp, Function
-        >>> f = Function('f')
-
         >>> solve(f(x) - x, f(x))
         [x]
         >>> solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
@@ -503,28 +470,26 @@ def solve(f, *symbols, **flags):
         >>> solve(x + exp(x)**2, exp(x), set=True)
         ([exp(x)], {(-sqrt(-x),), (sqrt(-x),)})
 
-        >>> from sympy import Indexed, IndexedBase, Tuple
+        >>> from sympy import Indexed, IndexedBase, Tuple, sqrt
         >>> A = IndexedBase('A')
         >>> eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
         >>> solve(eqs, eqs.atoms(Indexed))
         {A[1]: 1, A[2]: 2}
 
-        * To solve for a function within a derivative, use :func:`~.dsolve`.
-
-    To solve for a symbol implicitly, use implicit=True:
+        * To solve for a symbol implicitly, use implicit=True:
 
-        >>> solve(x + exp(x), x)
-        [-LambertW(1)]
-        >>> solve(x + exp(x), x, implicit=True)
-        [-exp(x)]
+            >>> solve(x + exp(x), x)
+            [-LambertW(1)]
+            >>> solve(x + exp(x), x, implicit=True)
+            [-exp(x)]
 
-    It is possible to solve for anything in an expression that can be
-    replaced with a symbol using :obj:`~sympy.core.basic.Basic.subs`:
+        * It is possible to solve for anything that can be targeted with
+          subs:
 
-        >>> solve(x + 2 + sqrt(3), x + 2)
-        [-sqrt(3)]
-        >>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2)
-        {y: -2 + sqrt(3), x + 2: -sqrt(3)}
+            >>> solve(x + 2 + sqrt(3), x + 2)
+            [-sqrt(3)]
+            >>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2)
+            {y: -2 + sqrt(3), x + 2: -sqrt(3)}
 
         * Nothing heroic is done in this implicit solving so you may end up
           with a symbol still in the solution:
@@ -535,7 +500,7 @@ def solve(f, *symbols, **flags):
             >>> solve(eqs, y*x, x)
             {x: -y - 4, x*y: -3*y - sqrt(3)}
 
-        * If you attempt to solve for a number, remember that the number
+        * If you attempt to solve for a number remember that the number
           you have obtained does not necessarily mean that the value is
           equivalent to the expression obtained:
 
@@ -546,22 +511,118 @@ def solve(f, *symbols, **flags):
             >>> [_.subs(z, -1) for _ in solve((x - y + 1).subs(-1, z), 1)]
             [-x + y]
 
+        * To solve for a function within a derivative, use ``dsolve``.
+
+    Single expression and more than one symbol:
+
+        * When there is a linear solution:
+
+            >>> solve(x - y**2, x, y)
+            [(y**2, y)]
+            >>> solve(x**2 - y, x, y)
+            [(x, x**2)]
+            >>> solve(x**2 - y, x, y, dict=True)
+            [{y: x**2}]
+
+        * If there is no linear solution, then the first successful
+          attempt for a nonlinear solution will be returned:
+
+            >>> solve(x**2 - y**2, x, y, dict=True)
+            [{x: -y}, {x: y}]
+            >>> solve(x**2 - y**2/exp(x), x, y, dict=True)
+            [{x: 2*LambertW(-y/2)}, {x: 2*LambertW(y/2)}]
+            >>> solve(x**2 - y**2/exp(x), y, x)
+            [(-x*sqrt(exp(x)), x), (x*sqrt(exp(x)), x)]
+
+        * When undetermined coefficients are identified:
+
+            This happens when it is possible to form a linear set of
+            equations in the variables provided from the coefficients
+            of the expressions in symbols not provided. A single
+            dictionary with specified values will be returned:
+
+            >>> eq = (a + b)*x - b + 2
+            >>> solve(eq, a, b)
+            {a: -2, b: 2}
+
+            The coefficient system solved was:
+
+            >>> list(eq.expand().as_coefficients_dict(x).values())
+            [a + b, 2 - b]
+
+            To obtain an algebraic solution in terms of ``a`` or ``b``
+            pass the equation in a list:
+
+            >>> solve([eq], a, b)
+            {a: b*(1 - x)/x - 2/x}
+
+    Iterable of one or more of the above:
+
+        * Involving relationals or bools:
+
+            >>> solve([x < 3, x - 2])
+            Eq(x, 2)
+            >>> solve([x > 3, x - 2])
+            False
+
+        * When the system is linear:
+
+            * With a solution:
+
+                >>> solve([x - 3], x)
+                {x: 3}
+                >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y)
+                {x: -3, y: 1}
+                >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y, z)
+                {x: -3, y: 1}
+                >>> solve((x + 5*y - 2, -3*x + 6*y - z), z, x, y)
+                {x: 2 - 5*y, z: 21*y - 6}
+
+            * Without a solution:
+
+                >>> solve([x + 3, x - 3])
+                []
+
+        * When the system is not linear:
+
+            >>> solve([x**2 + y -2, y**2 - 4], x, y, set=True)
+            ([x, y], {(-2, -2), (0, 2), (2, -2)})
+
+        * If no *symbols* are given, all free *symbols* will be selected and a
+          list of mappings returned:
+
+            >>> solve([x - 2, x**2 + y])
+            [{x: 2, y: -4}]
+            >>> solve([x - 2, x**2 + f(x)], {f(x), x})
+            [{x: 2, f(x): -4}]
+
+        * If any equation does not depend on the symbol(s) given, it will be
+          eliminated from the equation set and an answer may be given
+          implicitly in terms of variables that were not of interest:
+
+            >>> solve([x - y, y - 3], x)
+            {x: y}
+
     **Additional Examples**
 
     ``solve()`` with check=True (default) will run through the symbol tags to
     eliminate unwanted solutions. If no assumptions are included, all possible
     solutions will be returned:
 
+        >>> from sympy import Symbol, solve
         >>> x = Symbol("x")
         >>> solve(x**2 - 1)
         [-1, 1]
 
-    By setting the ``positive`` flag, only one solution will be returned:
+    By using the positive tag, only one solution will be returned:
 
         >>> pos = Symbol("pos", positive=True)
         >>> solve(pos**2 - 1)
         [1]
 
+    Assumptions are not checked when ``solve()`` input involves
+    relationals or bools.
+
     When the solutions are checked, those that make any denominator zero
     are automatically excluded. If you do not want to exclude such solutions,
     then use the check=False option:
@@ -570,15 +631,15 @@ def solve(f, *symbols, **flags):
         >>> solve(sin(x)/x)  # 0 is excluded
         [pi]
 
-    If ``check=False``, then a solution to the numerator being zero is found
-    but the value of $x = 0$ is a spurious solution since $\sin(x)/x$ has the well
-    known limit (without discontinuity) of 1 at $x = 0$:
+    If check=False, then a solution to the numerator being zero is found: x = 0.
+    In this case, this is a spurious solution since $\sin(x)/x$ has the well
+    known limit (without dicontinuity) of 1 at x = 0:
 
         >>> solve(sin(x)/x, check=False)
         [0, pi]
 
     In the following case, however, the limit exists and is equal to the
-    value of $x = 0$ that is excluded when check=True:
+    value of x = 0 that is excluded when check=True:
 
         >>> eq = x**2*(1/x - z**2/x)
         >>> solve(eq, x)
@@ -590,63 +651,6 @@ def solve(f, *symbols, **flags):
         >>> limit(eq, x, 0, '+')
         0
 
-    **Solving Relationships**
-
-    When one or more expressions passed to ``solve`` is a relational,
-    a relational result is returned (and the ``dict`` and ``set`` flags
-    are ignored):
-
-        >>> solve(x < 3)
-        (-oo < x) & (x < 3)
-        >>> solve([x < 3, x**2 > 4], x)
-        ((-oo < x) & (x < -2)) | ((2 < x) & (x < 3))
-        >>> solve([x + y - 3, x > 3], x)
-        (3 < x) & (x < oo) & Eq(x, 3 - y)
-
-    Although checking of assumptions on symbols in relationals
-    is not done, setting assumptions will affect how certain
-    relationals might automatically simplify:
-
-        >>> solve(x**2 > 4)
-        ((-oo < x) & (x < -2)) | ((2 < x) & (x < oo))
-
-        >>> r = Symbol('r', real=True)
-        >>> solve(r**2 > 4)
-        (2 < r) | (r < -2)
-
-    There is currently no algorithm in SymPy that allows you to use
-    relationships to resolve more than one variable. So the following
-    does not determine that ``q < 0`` (and trying to solve for ``r``
-    and ``q`` will raise an error):
-
-        >>> from sympy import symbols
-        >>> r, q = symbols('r, q', real=True)
-        >>> solve([r + q - 3, r > 3], r)
-        (3 < r) & Eq(r, 3 - q)
-
-    You can directly call the routine that ``solve`` calls
-    when it encounters a relational: :func:`~.reduce_inequalities`.
-    It treats Expr like Equality.
-
-        >>> from sympy import reduce_inequalities
-        >>> reduce_inequalities([x**2 - 4])
-        Eq(x, -2) | Eq(x, 2)
-
-    If each relationship contains only one symbol of interest,
-    the expressions can be processed for multiple symbols:
-
-        >>> reduce_inequalities([0 <= x  - 1, y < 3], [x, y])
-        (-oo < y) & (1 <= x) & (x < oo) & (y < 3)
-
-    But an error is raised if any relationship has more than one
-    symbol of interest:
-
-        >>> reduce_inequalities([0 <= x*y  - 1, y < 3], [x, y])
-        Traceback (most recent call last):
-        ...
-        NotImplementedError:
-        inequality has more than one symbol of interest.
-
     **Disabling High-Order Explicit Solutions**
 
     When solving polynomial expressions, you might not want explicit solutions


Test patch:
diff --git a/sympy/integrals/tests/test_risch.py b/sympy/integrals/tests/test_risch.py
index 51dd884a05..68be260e1f 100644
--- a/sympy/integrals/tests/test_risch.py
+++ b/sympy/integrals/tests/test_risch.py
@@ -374,6 +374,7 @@ def test_integrate_hyperexponential_returns_piecewise():
     # TODO: Add a test where two different parts of the extension use a
     # Piecewise, like y**x + z**x.
 
+
 def test_issue_13947():
     a, t, s = symbols('a t s')
     assert risch_integrate(2**(-pi)/(2**t + 1), t) == \
@@ -381,6 +382,7 @@ def test_issue_13947():
     assert risch_integrate(a**(t - s)/(a**t + 1), t) == \
         exp(-s*log(a))*log(a**t + 1)/log(a)
 
+
 def test_integrate_primitive():
     DE = DifferentialExtension(extension={'D': [Poly(1, x), Poly(1/x, t)],
         'Tfuncs': [log]})
@@ -738,3 +740,24 @@ def test_DifferentialExtension_printing():
     assert str(DE) == ("DifferentialExtension({fa=Poly(t1 + t0**2, t1, domain='ZZ[t0]'), "
         "fd=Poly(1, t1, domain='ZZ'), D=[Poly(1, x, domain='ZZ'), Poly(2*x*t0, t0, domain='ZZ[x]'), "
         "Poly(2*t0*x/(t0 + 1), t1, domain='ZZ(x,t0)')]})")
+
+
+def test_issue_23948():
+    f = (
+        ( (-2*x**5 + 28*x**4 - 144*x**3 + 324*x**2 - 270*x)*log(x)**2
+         +(-4*x**6 + 56*x**5 - 288*x**4 + 648*x**3 - 540*x**2)*log(x)
+         +(2*x**5 - 28*x**4 + 144*x**3 - 324*x**2 + 270*x)*exp(x)
+         +(2*x**5 - 28*x**4 + 144*x**3 - 324*x**2 + 270*x)*log(5)
+         -2*x**7 + 26*x**6 - 116*x**5 + 180*x**4 + 54*x**3 - 270*x**2
+        )*log(-log(x)**2 - 2*x*log(x) + exp(x) + log(5) - x**2 - x)**2
+       +( (4*x**5 - 44*x**4 + 168*x**3 - 216*x**2 - 108*x + 324)*log(x)
+         +(-2*x**5 + 24*x**4 - 108*x**3 + 216*x**2 - 162*x)*exp(x)
+         +4*x**6 - 42*x**5 + 144*x**4 - 108*x**3 - 324*x**2 + 486*x
+        )*log(-log(x)**2 - 2*x*log(x) + exp(x) + log(5) - x**2 - x)
+    )/(x*exp(x)**2*log(x)**2 + 2*x**2*exp(x)**2*log(x) - x*exp(x)**3
+       +(-x*log(5) + x**3 + x**2)*exp(x)**2)
+
+    F = ((x**4 - 12*x**3 + 54*x**2 - 108*x + 81)*exp(-2*x)
+        *log(-x**2 - 2*x*log(x) - x + exp(x) - log(x)**2 + log(5))**2)
+
+    assert risch_integrate(f, x) == F
diff --git a/sympy/physics/continuum_mechanics/tests/test_beam.py b/sympy/physics/continuum_mechanics/tests/test_beam.py
index 2c33fca5f9..f8b5ef202a 100644
--- a/sympy/physics/continuum_mechanics/tests/test_beam.py
+++ b/sympy/physics/continuum_mechanics/tests/test_beam.py
@@ -767,16 +767,3 @@ def test_max_deflection_Beam3D():
     p = sympify("-10 + 10*sqrt(10793)/43")
     q = sympify("(10 - 10*sqrt(10793)/43)**3/160 - 20/7 + (10 - 10*sqrt(10793)/43)**4/6400 + 20*sqrt(10793)/301 + 27*(10 - 10*sqrt(10793)/43)**2/560")
     assert b.max_deflection() == [(0, 0), (10, c), (p, q)]
-
-def test_torsion_Beam3D():
-    x = symbols('x')
-    b = Beam3D(20, 40, 21, 100, 25)
-    b.apply_moment_load(15, 5, -2, dir='x')
-    b.apply_moment_load(25, 10, -2, dir='x')
-    b.apply_moment_load(-5, 20, -2, dir='x')
-    b.solve_for_torsion()
-    assert b.angular_deflection().subs(x, 3) == sympify("1/40")
-    assert b.angular_deflection().subs(x, 9) == sympify("17/280")
-    assert b.angular_deflection().subs(x, 12) == sympify("53/840")
-    assert b.angular_deflection().subs(x, 17) == sympify("2/35")
-    assert b.angular_deflection().subs(x, 20) == sympify("3/56")
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 2f925c855b..02a81bffeb 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -2474,10 +2474,10 @@ def test_issue_20902():
 def test_issue_21034():
     a = symbols('a', real=True)
     system = [x - cosh(cos(4)), y - sinh(cos(a)), z - tanh(x)]
-    # constants inside hyperbolic functions should not be rewritten in terms of exp
     assert solve(system, x, y, z) == [(cosh(cos(4)), sinh(cos(a)), tanh(cosh(cos(4))))]
-    # but if the variable of interest is present in a hyperbolic function,
-    # then it should be rewritten in terms of exp and solved further
+    # constants inside hyperbolic functions should not be rewritten
+    # in terms of exp; rewriting should only happen with hyperbolics
+    # containing a symbol if interest
     newsystem = [(exp(x) - exp(-x)) - tanh(x)*(exp(x) + exp(-x)) + x - 5]
     assert solve(newsystem, x) == {x: 5}