While the curse of high dimensionality associated with the Fokker-Planck equation can be resolved by Monte Carlo solution to the underlying stochastic process, the distribution of discretization points for the resulting random walk-based method follows the evolution of the target distribution function, leading to the weight-spreading phenomenon in the $\delta f$ method, and therefore increase in statistical sampling error. In order to overcome this issue, we propose a Monte Carlo solution algorithm to the Fokker-Planck equation with fixed particles in phase space using the method of fundamental solutions. The numerical challenges introduced by the integral solution in resolving small time step sizes are addressed and dealt with by approximating Dirac delta functions with Gaussian kernels of finite variance around each particle. The numerical scheme is tested against the benchmark by studying the relaxation of non-equilibrium distribution functions following Brownian motion and the linearized Landau equation. For the relevant and finite range of time step sizes, the proposed convolution solution method provides a reasonable accuracy, yet the limit of its capability is investigated.