A numerical method, suitable for the simulation of the time evolution of quantum spin models of arbitrary lattice dimension, is presented. The method combines sampling of the Wigner function with evolution equations obtained from the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. Going to higher orders of the BBGKY hierarchy allows for a systematic refinement of the method. Quantum correlations are treated through both, the Wigner function sampling and the BBGKY evolution, bringing about highly accurate estimates of correlation functions. The method is particularly suitable for long-range interacting systems, and we demonstrate its power by comparing with exact results as well as other numerical methods. As an application we compute spin squeezing in a two-dimensional lattice with power-law interactions and a transverse field, which should be accessible in future ion trap experiments.