Abstract: In quasiclassical Keldysh theory, the Green function matrix g<cb><2021> is used to compute a variety of physical quantities in mesoscopic systems. However, solving the set of non-linear differential equations that provide g<cb><2021> becomes a challenging task when going to higher spatial dimensions than one. Such an extension is crucial in order to describe physical phenomena like charge/spin Hall effects and topological excitations like vortices and skyrmions, none of which can be captured in one-dimensional models. We here present a numerical finite element method which solves the 2D and 3D quasiclassical Usadel equation, without any linearisation, relevant for the diffusive regime. We show the application of this on two model systems with non-trivial geometries: (i) a bottlenecked Josephson junction with external flux and (ii) a nanodisk ferromagnet deposited on top of a superconductor. We demonstrate that it is possible to control externally not only the geometrical array in which superconducting vortices arrange themselves, but also to cause coalescence and thus tune the number of vortices. The finite element method presented herein could pave the way for gaining insight in physical phenomena which so far have remained largely unexplored due to the complexity of solving the full quasiclassical equations in higher dimensions.