We describe a general Godunov-type splitting for numerical simulations of the Fisher-Kolmogorov-Petrovski-Piskunov growth and diffusion equation in two spatial dimensions with Neumann boundary conditions. In particular, the method is appropriate for modeling population growth and dispersal on a terrestrial map. The procedure is semi-implicit, hence quite stable, and approximately O(x^2)+O(t^2) accurate, excluding boundary condition complications. It also has low memory requirements and shows good performance. Our principal application of this solver is global human dispersal in the late Pleistocene, modeled via growth and diffusion over geographical maps of changing paleovegetation and paleoclimate.