Many times in quantum physics we face a problem of calculating expectation values of normal product operators. Assume we have a two level system with creation and annihilation operators $\hat{a}$ and $\hat{b}$. We would like to find the time evolution of expectation values of the form
$$\langle (\hat{a}^{\dagger})^{k}(\hat{b})^{l}\rangle$$ $$\langle (\hat{b}^{\dagger})^{k}(\hat{a})^{l}\rangle$$
We can do that using generating function method. For example, if we take Hamiltonian $H = S_{z}^2 = \frac{1}{4}(\hat{a}^{\dagger}\hat{a} - \hat{b}^{\dagger}\hat{b})^2$, equation for density matrix coefficients expanded in the Fock state basis become: $$\hat{\rho}(t) = \sum\limits_{k_1,k_2,n_1,n_2} \rho^{k_1,k_2}_{n_1,n_2}(t) |k_1,k_2\rangle\langle n_1,n_2|.$$ using von Neumann equation we arrive at: $$i\dot{\rho}^{k_1,k_2}_{n_1,n_2}(t) = \frac{1}{4}\rho^{k_1,k_2}_{n_1,n_2}(t)\left[ (k_1-k_2)^2 - (n_1-n_2)^2\right]$$ Of course we can solve this equation easily, but the idea is different. Let's say we wan to calculate expectation value: $$\langle \hat{a}^{\dagger}\hat{b}\rangle = \sum\limits_{k_1,k_2}\rho^{k_1-1,k_2+1}_{k_1,k_2}\sqrt{k_1(k_2+1)}$$ We introduce generating function of the form: $$G(x,y,t) = \sum\limits_{k_1,k_2}x^{k_1}y^{k_2}\rho^{k_1-1,k_2+1}_{k_1,k_2}(t)\sqrt{k_1(k_2+1)}$$ such that $G(1,1,t) = \langle \hat{a}^{\dagger}\hat{b}\rangle$. Using recurrence relation we can write $$i\partial_{t}G(x,y,t) = (1-x\partial_{x} + y\partial_{y})G(x,y,t)$$ We can solve the differential equation and get the result.
This was a simple example just to show the idea. What about more complicated Hamiltonians like: $H = S_{x} = \frac{1}{2}(\hat{a}^{\dagger}\hat{b} + \hat{b}^{\dagger}\hat{a})$ or $H = S_{z}S_{y} + S_{y}S_{z}$, where $S_y = \frac{1}{2i}(\hat{a}^{\dagger}\hat{b} - \hat{b}^{\dagger}\hat{a})$.
Is it a general trick, or it works only for diagonal Hamiltonians? Maybe some of you know other methods that can be used for calculating expectation values of operator products?