A numerical method is developed for approximating the solution to the Vlasov-Poisson-Fokker-Planck system in two spatial dimensions. The method generalizes the approximation for the system in one dimension given in [S. Wollman, E. Ozizmir, Numerical approximation of the Vlasov-Poisson-Fokker-Planck system in one dimension, J. Comput. Phys. 202 (2005) 602-644]. The numerical procedure is based on a change of variables that puts the convection-diffusion equation into a form so that finite difference methods for parabolic type partial differential equations can be applied. The computational cycle combines a type of deterministic particle method with a periodic interpolation of the solution along particle trajectories onto a fixed grid.computational work is done to demonstrate the accuracy and effectiveness of the approximation method. Parts of the numerical procedure are adapted to run on a parallel computer.