title: Solving system of polynomial matrix equations in sympy
body: <p>I would like to use sympy to solve for a vector that satisfies a set of quadratic and linear equations, such as <code>x.T @ A @ x = x.T @ B @ x = 1</code>, and <code>x.T @ C = 0</code>. Here is my attempt:</p>
  <pre><code>import numpy as np
  from sympy.polys.polymatrix import PolyMatrix
  from sympy import symbols
  from sympy.solvers import solve_poly_system
  
  x = PolyMatrix([symbols(f'x{k}') for k in range(2)])
  A = np.eye(2)
  B = np.array([1, 0]).reshape([2, 1])
  one = PolyMatrix([1])
  
  solve_poly_system([x.T @ A @ x - one, x.T @ B], [symbols(f'x{k}') for k in range(2)])
  </code></pre>
  <p>However, this generates the following error:</p>
  <pre><code>PolynomialError: non-commutative expressions are not supported
  </code></pre>
  <p>Any suggestions?</p>