A method is presented to solve the Bogoliubov-de Gennes equations with arbitrary distributions of vortices. The real-space Green's function approach based on Chebyshev polynomials is complemented by a gauge transformation which allows one to treat finite as well as infinite, ordered as well as disordered vortex configurations. This tool gives unprecedented access to vortex lattices at very low magnetic fields and glassy phases. After describing in detail the method and its implementation, we use it to address a series of problems related to d-wave superconductivity on the square lattice. We first study the continuity of the vortex-core energy spectrum and its evolution from the quantum regime to the semiclassical limit; we investigate the effect of the band structure on the vortex by following the self-consistent solution through a Lifshitz transition; we then study the evolution from the vortex lattice to the isolated-vortex limit with decreasing field and show that a new emerging length scale controls this transition; finally, we perform a statistical study of the vortex-core local density of states in the presence of positional disorder in the vortex lattice. The calculations reveal a number of qualitative differences between the properties of vortices in the quantum and semiclassical regimes.