matrix - How to force Series[] to compute expansions by considering non commutative multiplication?

I wish to compute the Taylor series expansion of the following iteration method $x_{k+1}=x_k-f'(x_k)^{-1}f(x_k)$ up to four terms of error ($e_k=x_k-\alpha$). When this is a scalar iteration, I very simply write the following

ClearAll["Global`*"]
f[e_] := df*(e^1 + c2 e^2 + c3 e^3 + c4 e^4 + c5 e^5);
fe = f[e];
f1e = f'[e];

Series[f1e^-1, {e, 0, 4}]*df // Simplify
u = e - Series[fe/f1e, {e, 0, 4}] // Simplify

and obtain correct results $$f(x_k)^{-1}=[1-2 \text{c2} e_k+ \left(4 \text{c2}^2-3 \text{c3}\right)e_k^2$$ $$-4 (2 c2^3 - 3 c2 c3 + c4) e_k^3 + (16 c2^4 - 36 c2^2 c3 + 9 c3^2 + 16 c2 c4 - 5 c5) e_k^4]f'({\alpha})^{-1}.$$ But, how to do this for the multi-dimensional case. That is, when the coefficients $c2,c3,c4$, and even $df=f'({\alpha})^{-1}$ are all matrices (note that e.g. $c2=\frac{1}{2!}f'(\alpha)^{-1}f^{(2)}(\alpha)$). By hand, I obtain the following correct results: $$ f(x_k)^{-1}=\left(I-2c2e_k+(4c2^2-3c3)e_k^2+(6c3c2+6c2c3-8c2^3-4c4)e_k^3\\ + (8c4c2+9c3^2+8c2c4-5c5-12c3c2^2-12c2c3c2-12c2^2c3+16c2^4)e_k^4\right)f'(\alpha)^{-1} $$ and \begin{equation} e_{k+1}=x_{k+1}-\alpha=-c2e_k^2+(2c2^2-2c3)e_k^3+(-4c2^3+4c2c3+3c3c2-3c4)e_k^4+\mathcal{O}(e_k^5). \end{equation}

So, my question is how to handle this non-commutative multiplication inside Series[] (e.g. $c2c3$ is not equal to $c3c2$ in this case)? I also tried to apply $**$, but I failed. Any suggestions or tricks to solve this problem will be appreciated fully.

I also saw the following posts in MathematicaStackExchange, but they were not useful for this problem, How to make noncommutative multiplication agree with commutative multiplication and Noncommutative multiply- expand expression.