I have been working on symbolic computations of bosonic algebra for decades now, see in particular posts on Blog:Notes, especially NormalOrder[] and examples of its use. Here I actualize some observations on the Mathematica piece of code (BosonNormalOrder) which I use to check such results.
One useful thing is to have a typeset version of the results. One can try to do that from within the notebook (a-Ctrl+^-esc-\dagger-esc, etc. for $a^\dagger$) but for complicated expressions, this is a poor approach as Mathematica is quite unstable and will revert (with no good undo facility) to some pseudocode of his. Therefore it is better to typeset it in $\mathrm{\LaTeX}$ and then convert that into some Mathematica form, e.g.,
ToExpression["aa^{\\dagger k}a^la^\\dagger=a^{\\dagger(k+1)}a^{l+1}+(k+l+1)a^{\\dagger k}a^l+kla^{\\dagger(k-1)}a^{l-1}", TeXForm, HoldForm]
where the bit between quotes is copy-pasted from Emacs (in which case Mathematica kindly offers to do the double escaping of slashes itself). This produces something that can be put in a text explaining what is going on:
And one has the main layout of how the test works: for algebraic values $k$ and $l$ (integers, possibly zero) but which are left undetermined, I haven't, yet, written it in a way that the calculation proceeds symbolically. It can definitely be done. For now, however, I worry more on making sure that all expressions I use are correct than automatizing everything. So I just check the sought result to my result, calculated by hand, and written in the form of a "wseq" (weighted sequence). If that matches for some particular cases, chances are, there are no mistakes. The two corresponding parts are the two grey lines.