Repository: sympy/sympy

Issue:
`solve_undetermined_coeffs` does not respect `dict=True`.
This is a regression in master since SymPy 1.10.1 (due to #23877).

The `solve_determined_coeffs` function no longer respects the `dict=True` option.

With 1.10.1:
```python
In [1]: a, b, c = symbols('a, b, c')

In [2]: solve_undetermined_coeffs(a*x + b-1, [a, b], [x])
Out[2]: {a: 0, b: 1}

In [3]: solve_undetermined_coeffs(a*x + b-1, [a, b], [x], dict=True)
Out[3]: [{a: 0, b: 1}]
```
With master:
```python
In [1]: a, b, c = symbols('a, b, c')

In [2]: solve_undetermined_coeffs(a*x + b-1, [a, b], [x])
Out[2]: {a: 0, b: 1}

In [3]: solve_undetermined_coeffs(a*x + b-1, [a, b], [x], dict=True)
Out[3]: {a: 0, b: 1}   <---- should be a list of dicts
```

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).



* No release notes entry will be added for this pull request.

<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. -->


    #### Brief description of what is fixed or changed


    #### 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 -->
    NO ENTRY
    <!-- END RELEASE NOTES -->


</details><p>

smichr: per #23927, consider returning empty {} when there is no solution -- will this work with `recurr` expectations which depend on getting a [] for something like `solve_undetermined_coeffs(0, [a,b,c], x)`?
smichr: check
```
File "/home/runner/work/sympy/sympy/sympy/solvers/tests/test_solvers.py", line 172, in test_solve_args
    assert solve(a*x + b**2/(x + 4) - 3*x - 4/x, a, b, dict=True) == \
AssertionError
________________________________________________________________________________
_____________ sympy/solvers/tests/test_solvers.py:test_issue_5901 ______________
Traceback (most recent call last):
  File "/home/runner/work/sympy/sympy/sympy/solvers/tests/test_solvers.py", line 1517, in test_issue_5901
    assert solve(-f(a)**2*g(a)**2 + f(a)**2*h(a)**2 + g(a).diff(a),
AssertionError: ([h(a), g(a)], {(h(a), g(a))})
________________________________________________________________________________
_____________ sympy/solvers/tests/test_solvers.py:test_issue_1[202](https://github.com/sympy/sympy/runs/7823958140?check_suite_focus=true#step:6:203)4 _____________
Traceback (most recent call last):
  File "/home/runner/work/sympy/sympy/sympy/solvers/tests/test_solvers.py", line 2404, in test_issue_12024
    assert solve(Piecewise((0.0, x < 0.1), (x, x >= 0.1)) - y) == \
AssertionError
```
sympy-bot: ### 🟠

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at https://github.com/sympy/sympy-bot/issues/75.

The following commits **add new files**:
* 7fb44844f74e3a20fd1beb1bb0639674ade2987b:
  - `doc/src/explanation/solve_output.rst`

If these files were added/deleted on purpose, you can ignore this message.

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
     [8ce582b4]       [f0154339]
-     1.20±0.01ms          763±5μs     0.64  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-     3.47±0.01ms      1.40±0.01ms     0.40  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     6.94±0.08ms      2.06±0.02ms     0.30  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).
smichr: - [ ] rebase on master to remove already committed lines
sympy-bot: ### 🟠

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at https://github.com/sympy/sympy-bot/issues/75.

The following commits **add new files**:
* 7fb44844f74e3a20fd1beb1bb0639674ade2987b:
  - `doc/src/explanation/solve_output.rst`

If these files were added/deleted on purpose, you can ignore this message.

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
     [8ce582b4]       [0f809ebc]
-         991±2μs          630±2μs     0.64  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-     2.85±0.01ms         1.16±0ms     0.41  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.69±0.02ms         1.71±0ms     0.30  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).
smichr: Any changes you would suggest, @oscarbenjamin ?
smichr: When the tests pass I will squash to a single commit and merge if there is no objection. This function is now just a front end to calling `solve` with flags after extracting the coefficient system.

Patch:
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/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..6b410e7874 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -1275,22 +1275,19 @@ def _has_piecewise(e):
 
 
 def _solve_undetermined(g, symbols, flags):
-    # helper to find independent variable for undetermined coefficients
-    # call and to return a list with a single dictionary else None
-    #                                              a*x + 2*x + b - c
-    symset = set(symbols)                        # {a, b}
-    free = g.free_symbols                        # {a, b, c, x}
-    ex = free - symset                           # {c, x}
-    if len(ex) != 1:
-        ind, dep = g.as_independent(*symbols)    # (2*x - c, a*x + b)
-        ex = ind.free_symbols & dep.free_symbols # {x, c} & {a, x, b} -> {x}
-    if len(ex) != 1:                             # e.g. (a + b)*x + b - c -> {c}, {a, b, x}
-        ex = dep.free_symbols - symset           # {x} = {a, b, x}-{a, b}
-    if len(ex) == 1:
-        ex = ex.pop()
-        sol = solve_undetermined_coeffs(g, symbols, ex, **flags)
-        if type(sol) is dict:
-            return [sol]
+    """solve helper to return a list with one dict (solution) else None
+
+    A direct call to solve_undetermined_coeffs is more flexible and
+    can return both multiple solutions and handle more than one independent
+    variable. Here, we have to be more cautious to keep from solving
+    something that does not look like an undetermined coeffs system --
+    to minimize the surprise factor since singularities that cancel are not
+    prohibited in solve_undetermined_coeffs.
+    """
+    if g.free_symbols - set(symbols):
+        sol = solve_undetermined_coeffs(g, symbols, **dict(flags, dict=True, set=None))
+        if len(sol) == 1:
+            return sol
 
 
 def _solve(f, *symbols, **flags):
@@ -2305,60 +2302,136 @@ def solve_linear_system(system, *symbols, **flags):
     return sol
 
 
-def solve_undetermined_coeffs(equ, coeffs, sym, **flags):
+def solve_undetermined_coeffs(equ, coeffs, *syms, **flags):
     r"""
-    Solve equation of a type $p(x; a_1, \ldots, a_k) = q(x)$ where both
-    $p$ and $q$ are univariate polynomials that depend on $k$ parameters.
+    Solve a system of equations in $k$ parameters that is formed by
+    matching coefficients in variables ``coeffs`` that are on
+    factors dependent on the remaining variables (or those given
+    explicitly by ``syms``.
 
     Explanation
     ===========
 
     The result of this function is a dictionary with symbolic values of those
-    parameters with respect to coefficients in $q$, [] if there is no
-    solution, else None if the system was not recognized.
-
-    This function accepts both equations class instances and ordinary
-    SymPy expressions. Specification of parameters and variables is
-    obligatory for efficiency and simplicity reasons.
+    parameters with respect to coefficients in $q$ -- empty if there
+    is no solution or coefficients do not appear in the equation -- else
+    None (if the system was not recognized). If there is more than one
+    solution, the solutions are passed as a list. The output can be modified using
+    the same semantics as for `solve` since the flags that are passed are sent
+    directly to `solve` so, for example the flag ``dict=True`` will always return a list
+    of solutions as dictionaries.
+
+    This function accepts both Equality and Expr class instances.
+    The solving process is most efficient when symbols are specified
+    in addition to parameters to be determined,  but an attempt to
+    determine them (if absent) will be made. If an expected solution is not
+    obtained (and symbols were not specified) try specifying them.
 
     Examples
     ========
 
     >>> from sympy import Eq, solve_undetermined_coeffs
-    >>> from sympy.abc import a, b, c, x
+    >>> from sympy.abc import a, b, c, h, p, k, x, y
 
-    >>> solve_undetermined_coeffs(Eq(2*a*x + a+b, x), [a, b], x)
+    >>> solve_undetermined_coeffs(Eq(a*x + a + b, x/2), [a, b], x)
     {a: 1/2, b: -1/2}
+    >>> solve_undetermined_coeffs(a - 2, [a])
+    {a: 2}
 
-    >>> solve_undetermined_coeffs(Eq(a*c*x + a+b, x), [a, b], x)
-    {a: 1/c, b: -1/c}
+    The equation can be nonlinear in the symbols:
 
-    >>> solve_undetermined_coeffs(0, [a], x)
-    []
+    >>> X, Y, Z = y, x**y, y*x**y
+    >>> eq = a*X + b*Y + c*Z - X - 2*Y - 3*Z
+    >>> coeffs = a, b, c
+    >>> syms = x, y
+    >>> solve_undetermined_coeffs(eq, coeffs, syms)
+    {a: 1, b: 2, c: 3}
 
-    >>> solve_undetermined_coeffs(a**2*x + 2*x + b + 3, [a, b], x) is None
-    True
+    And the system can be nonlinear in coefficients, too, but if
+    there is only a single solution, it will be returned as a
+    dictionary:
+
+    >>> eq = a*x**2 + b*x + c - ((x - h)**2 + 4*p*k)/4/p
+    >>> solve_undetermined_coeffs(eq, (h, p, k), x)
+    {h: -b/(2*a), k: (4*a*c - b**2)/(4*a), p: 1/(4*a)}
 
+    Multiple solutions are always returned in a list:
+
+    >>> solve_undetermined_coeffs(a**2*x + b - x, [a, b], x)
+    [{a: -1, b: 0}, {a: 1, b: 0}]
+
+    Using flag ``dict=True`` (in keeping with semantics in :func:`~.solve`)
+    will force the result to always be a list with any solutions
+    as elements in that list.
+
+    >>> solve_undetermined_coeffs(a*x - 2*x, [a], dict=True)
+    [{a: 2}]
     """
+    if not (coeffs and all(i.is_Symbol for i in coeffs)):
+        raise ValueError('must provide symbols for coeffs')
+
     if isinstance(equ, Eq):
-        # got equation, so move all the
-        # terms to the left hand side
-        equ = equ.lhs - equ.rhs
+        eq = equ.lhs - equ.rhs
+    else:
+        eq = equ
+
+    ceq = cancel(eq)
+    xeq = _mexpand(ceq.as_numer_denom()[0], recursive=True)
+
+    free = xeq.free_symbols
+    coeffs = free & set(coeffs)
+    if not coeffs:
+        return ([], {}) if flags.get('set', None) else []  # solve(0, x) -> []
+
+    if not syms:
+        # e.g. A*exp(x) + B - (exp(x) + y) separated into parts that
+        # don't/do depend on coeffs gives
+        # -(exp(x) + y), A*exp(x) + B
+        # then see what symbols are common to both
+        # {x} = {x, A, B} - {x, y}
+        ind, dep = xeq.as_independent(*coeffs, as_Add=True)
+        dfree = dep.free_symbols
+        syms = dfree & ind.free_symbols
+        if not syms:
+            # but if the system looks like (a + b)*x + b - c
+            # then {} = {a, b, x} - c
+            # so calculate {x} = {a, b, x} - {a, b}
+            syms = dfree - set(coeffs)
+        if not syms:
+            syms = [Dummy()]
+    else:
+        if len(syms) == 1 and iterable(syms[0]):
+            syms = syms[0]
+        e, s, _ = recast_to_symbols([xeq], syms)
+        xeq = e[0]
+        syms = s
 
-    equ = cancel(equ).as_numer_denom()[0]
+    # find the functional forms in which symbols appear
 
-    system = list(collect(equ.expand(), sym, evaluate=False).values())
+    gens = set(xeq.as_coefficients_dict(*syms).keys()) - {1}
+    cset = set(coeffs)
+    if any(g.has_xfree(cset) for g in gens):
+        return  # a generator contained a coefficient symbol
 
-    if not any(equ.has(sym) for equ in system):
-        # consecutive powers in the input expressions have
-        # been successfully collected, so solve remaining
-        # system using Gaussian elimination algorithm
-        _flags = dict(flags, set=None, dict=True)
-        sol = solve(system, *coeffs, **_flags)
-        if len(sol) == 1:
-            return sol[0]
-        elif not sol:
-            return []
+    # make sure we are working with symbols for generators
+
+    e, gens, _ = recast_to_symbols([xeq], list(gens))
+    xeq = e[0]
+
+    # collect coefficients in front of generators
+
+    system = list(collect(xeq, gens, evaluate=False).values())
+
+    # get a solution
+
+    soln = solve(system, coeffs, **flags)
+
+    # unpack unless told otherwise if length is 1
+
+    settings = flags.get('dict', None) or flags.get('set', None)
+    if type(soln) is dict or settings or len(soln) != 1:
+        return soln
+    return soln[0]
 
 
 def solve_linear_system_LU(matrix, syms):
@@ -2533,12 +2606,12 @@ def _tsolve(eq, sym, **flags):
     Examples
     ========
 
-    >>> from sympy import log
+    >>> from sympy import log, ordered
     >>> from sympy.solvers.solvers import _tsolve as tsolve
     >>> from sympy.abc import x
 
-    >>> set(tsolve(3**(2*x + 5) - 4, x))
-    {(-5*log(3)/2 + log(2) + I*pi)/log(3), -5/2 + log(2)/log(3)}
+    >>> list(ordered(tsolve(3**(2*x + 5) - 4, x)))
+    [-5/2 + log(2)/log(3), (-5*log(3)/2 + log(2) + I*pi)/log(3)]
 
     >>> tsolve(log(x) + 2*x, x)
     [LambertW(2)/2]


Test patch:
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..7bea33398a 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -736,15 +736,43 @@ def test_solve_linear():
 
 
 def test_solve_undetermined_coeffs():
-    assert solve_undetermined_coeffs(a*x**2 + b*x**2 + b*x + 2*c*x + c + 1, [a, b, c], x) == \
-        {a: -2, b: 2, c: -1}
+    assert solve_undetermined_coeffs(
+        a*x**2 + b*x**2 + b*x + 2*c*x + c + 1, [a, b, c], x
+        ) == {a: -2, b: 2, c: -1}
     # Test that rational functions work
-    assert solve_undetermined_coeffs(a/x + b/(x + 1) - (2*x + 1)/(x**2 + x), [a, b], x) == \
-        {a: 1, b: 1}
+    assert solve_undetermined_coeffs(a/x + b/(x + 1)
+        - (2*x + 1)/(x**2 + x), [a, b], x) == {a: 1, b: 1}
     # Test cancellation in rational functions
-    assert solve_undetermined_coeffs(((c + 1)*a*x**2 + (c + 1)*b*x**2 +
-    (c + 1)*b*x + (c + 1)*2*c*x + (c + 1)**2)/(c + 1), [a, b, c], x) == \
+    assert solve_undetermined_coeffs(
+        ((c + 1)*a*x**2 + (c + 1)*b*x**2 +
+        (c + 1)*b*x + (c + 1)*2*c*x + (c + 1)**2)/(c + 1),
+        [a, b, c], x) == \
         {a: -2, b: 2, c: -1}
+    # multivariate
+    X, Y, Z = y, x**y, y*x**y
+    eq = a*X + b*Y + c*Z - X - 2*Y - 3*Z
+    coeffs = a, b, c
+    syms = x, y
+    assert solve_undetermined_coeffs(eq, coeffs) == {
+        a: 1, b: 2, c: 3}
+    assert solve_undetermined_coeffs(eq, coeffs, syms) == {
+        a: 1, b: 2, c: 3}
+    assert solve_undetermined_coeffs(eq, coeffs, *syms) == {
+        a: 1, b: 2, c: 3}
+    # check output format
+    assert solve_undetermined_coeffs(a*x + a - 2, [a]) == []
+    assert solve_undetermined_coeffs(a**2*x - 4*x, [a]) == [
+        {a: -2}, {a: 2}]
+    assert solve_undetermined_coeffs(0, [a]) == []
+    assert solve_undetermined_coeffs(0, [a], dict=True) == []
+    assert solve_undetermined_coeffs(0, [a], set=True) == ([], {})
+    assert solve_undetermined_coeffs(1, [a]) == []
+    abeq = a*x - 2*x + b - 3
+    s = {b, a}
+    assert solve_undetermined_coeffs(abeq, s, x) == {a: 2, b: 3}
+    assert solve_undetermined_coeffs(abeq, s, x, set=True) == ([a, b], {(2, 3)})
+    assert solve_undetermined_coeffs(sin(a*x) - sin(2*x), (a,)) is None
+    assert solve_undetermined_coeffs(a*x + b*x - 2*x, (a, b)) == {a: 2 - b}
 
 
 def test_solve_inequalities():
@@ -2543,3 +2571,16 @@ def test_issue_10169():
         e: Rational(-40,129) - 5*2**Rational(1,4)/1032 + 5*2**Rational(3,4)/129,
         k: -10*sqrt(2)/129 + 5*2**Rational(1,4)/258 + 20*2**Rational(3,4)/129
     }
+
+
+def test_solve_undetermined_coeffs_issue_23927():
+    A, B, r, phi = symbols('A, B, r, phi')
+    eq = Eq(A*sin(t) + B*cos(t), r*sin(t - phi)).rewrite(Add).expand(trig=True)
+    soln = solve_undetermined_coeffs(eq, (r, phi), t)
+    assert soln == [{
+        phi: 2*atan((A - sqrt(A**2 + B**2))/B),
+        r: (-A**2 + A*sqrt(A**2 + B**2) - B**2)/(A - sqrt(A**2 + B**2))
+        }, {
+        phi: 2*atan((A + sqrt(A**2 + B**2))/B),
+        r: (A**2 + A*sqrt(A**2 + B**2) + B**2)/(A + sqrt(A**2 + B**2))/-1
+        }]