Given an arbitrary set of (numerical) square complex matrices $\mathcal{A}=\{A_1,A_2,\cdots,A_m\}$, I am interested in computing the real matrix Lie algebra generated by $\mathcal{A}$, call it $\mathcal{L_\mathcal{A}}$. That is, I would like a basis for $$ \mathcal{L_\mathcal{A}} = \mathbb{span_R}\{B:B\in\cup_{k=1}^{\infty}\mathcal{C}_k\} $$ where $\mathcal{C}_k$ is defined recursively as $\mathcal{C_1}=\mathcal{A}$, and $\mathcal{C_{k+1}}=\{[X,Y]:X,Y\in\cup_{j=1}^k\mathcal{C_j}\}$ for $k\geq 1$.
This calculation comes up in (quantum) control theory.
Currently I am using a method found here which searches only through repeated Lie brackets (i.e. ones of the form $[A_{j_1},[A_{j_2},[A_{j_3},\cdots[A_{j_{n-1}},A_{j_n}]\cdots]]]$), and is guaranteed to terminate. However I'm interested to know if there are any other (faster) methods. Perhaps using P. Hall bases? Perhaps a recursive algorithm? My default language at the moment is Matlab.