How to find noncommutative product of two expression

Question

I have the following function $$ B(n)=x \alpha ^n-y \alpha ^n q^n + \epsilon \left(x \alpha ^{n+1}-y \alpha ^{n+1} q^{n+1}\right). $$ Here, I want to find the product $B(n)^2$ or $B(n)*B(n+1)$ where the variables $x$ and $y$ are non-commutative i.e. $xy\neq yx$. Is there any way to achieve this in Mathematica ?

For example, when I try to find the product in Mathematica, I find that

$$ B(n)^2 = -2 x y \epsilon ^2 \alpha ^{2 n+2} q^{n+1}-2 x y \epsilon \alpha ^{2 n+1} q^{n+1}+y^2 \epsilon ^2 \alpha ^{2 n+2} q^{2 n+2}+2 y^2 \epsilon \alpha ^{2 n+1} q^{2 n+1}+y^2 \alpha ^{2 n} q^{2 n}-2 x y \epsilon \alpha ^{2 n+1} q^n-2 x y \alpha ^{2 n} q^n+x^2 \epsilon ^2 \alpha ^{2 n+2}+2 x^2 \epsilon \alpha ^{2 n+1}+x^2 \alpha ^{2 n} $$ and $$ B(n)*B(n+1) = -x y \epsilon ^2 \alpha ^{2 n+3} q^{n+1}-x y \epsilon ^2 \alpha ^{2 n+3} q^{n+2}-2 x y \epsilon \alpha ^{2 n+2} q^{n+1}-x y \epsilon \alpha ^{2 n+2} q^{n+2}-x y \alpha ^{2 n+1} q^{n+1}+y^2 \epsilon ^2 \alpha ^{2 n+3} q^{2 n+3}+2 y^2 \epsilon \alpha ^{2 n+2} q^{2 n+2}+y^2 \alpha ^{2 n+1} q^{2 n+1}-x y \epsilon \alpha ^{2 n+2} q^n-x y \alpha ^{2 n+1} q^n+x^2 \epsilon ^2 \alpha ^{2 n+3}+2 x^2 \epsilon \alpha ^{2 n+2}+x^2 \alpha ^{2 n+1}. $$ But it must be that $xy \neq yx$.

Answer 1

In the light of Mathematica's NonCommutativeMultiply, I have solved the problem as follows:

B[n_] := (a ** x + b ** y) + (b ** x + a ** y ** q^n);

ExpandNCM[(h : NonCommutativeMultiply)[a___, b_Plus, c___]] := 
 Distribute[h[a, b, c], Plus, h, Plus, ExpandNCM[h[##]] &]

ExpandNCM[(h : NonCommutativeMultiply)[a___, b_Times, c___]] := 
 Most[b] ExpandNCM[h[a, Last[b], c]]

ExpandNCM[a_] := ExpandAll[a]

B[n]**B[n + 1] // ExpandNCM // ExpandAll

= a ** x ** a ** x + a ** x ** b ** x + a ** x ** b ** y + 
  b ** x ** a ** x + b ** x ** b ** x + b ** x ** b ** y + 
  b ** y ** a ** x + b ** y ** b ** x + b ** y ** b ** y + 
  a ** x ** a ** y ** q^(1 + n) + a ** y ** q^n ** a ** x + 
  a ** y ** q^n ** b ** x + a ** y ** q^n ** b ** y + 
  b ** x ** a ** y ** q^(1 + n) + b ** y ** a ** y ** q^(1 + n) + 
  a ** y ** q^n ** a ** y ** q^(1 + n)

The only problem is that you must simplify some expressions with manually.

Thank you @CATrevillian.