Starting point for a semi-analytical solution of the Fokker-Planck forward equation for the Heston model is the exact sampling algorithm of Broadie and Kaya  (for the notation please see )
The probability distribution function can be described as
and is given by a noncentral chi squared distribution. The distribution can be calculated using the exact simulation algorithm. In this algorithm the variable is given as a function of two random variables and .
The distribution of can now be derived using the general transformation theorem for random variables: Let X be a random variable with probability density function f. The transformed random variable has the probability density function
First step is now to rewrite the exact simulation method in terms of the two random variables
The simulation scheme then becomes
or in terms of the transformed random variable
Let be the density function of
and follows by definition a normal distribution. The joint probability density function of ) is then
This yields to the semi-analytical formula for the solution of the Fokker-Planck equation because by definition is the distribution density function of , which is given by
The integration over can be carried out using e.g. the Simpson integral rule together with the Cornish-Fisher expansion, which gives an upper bound for the truncation of the upper limit of the integration.
The contour plots below show the probability density function of the Heston model for some example parametrisations.
The example code is available here and depends on the upcoming QuantLib version 1.4.
 M. Broadie, Ö. Kaya, Exact Simulation of Stochastic Volatility and other Affine Jump Diffusion Processes
 K. Spanderen, Probability Distribution of the Heston Model, Part I
 R.U. Seydel, Tools for Computational Finance, pp 86