I would like to visualize $W(X,P)=W(\alpha,\alpha^*)$ from a given density matrix $\hat{\rho}$ that has been obtained before e.g. from the master equation. I am especially interested in density matrices describing a single bosonic mode with Kerr nonlinearity and expressed in Fock states up to some cutoff $N_{\text{cutoff}}\sim 50$. The states are expected to be quite similar to a gaussian state (as a matter of fact, when it is known that the state is truly Gaussian, the Wigner function can readily be obtained from correlation functions, but I explicitly want to see how far this gaussian approximation is from the true state.)
My initial take at calculating the wigner function is using the Characteristic function: $$\chi_W(\lambda,\lambda^*)=Tr\left\{\rho \exp(\lambda a^\dagger-\lambda^*a) \right\}\\ W(\alpha,\alpha^*)=\frac{1}{\pi^2}\int d^2\lambda \exp(-\lambda\alpha^*+\lambda^*\alpha)\,\chi_W(\lambda,\lambda^*) $$ Unfortunately, for let's say a 100x100 grid in the $\alpha$- and $\lambda$-planes the integral sums over 10^8 terms, which takes ages (a day) to evaluate to get a result that is of a too low quality.
I found here a program which uses an alternative approach based on Laguerre polynomials (I haven't found where their algorithm stems from unfortunately). This is already significantly faster (~30 minutes) and gives a meaningful, smooth 'mountain' in phase space corresponding to the state of the system, but aside from that also a giant noisy ring around the origin of fluctuations (I thought for a while this could be physical interference, but it is too irregular and doesn't disappear when feeding a $\hat{\rho}$ that is explicitly gaussian thus classical).
The last source I found is this github page, but it is not entirely transparant to me (partly because I'm unfamiliar with the language). It does tell me that FFT is also an option, though I don't know how to do that for an integral over the whole complex plane.
I prefer to work with MATLAB, but Mathematica is also an option.
