An extension of quasiclassical Keldysh-Usadel theory to higher spatial dimensions than one is crucial in order to describe physical phenomena like charge/spin Hall effects and topological excitations like vortices and skyrmions, none of which are captured in one-dimensional models. We here present a numerical finite element method which solves the non-linearized 2D and 3D quasiclassical Usadel equation relevant for the diffusive regime. We show the application of this on three model systems with non-trivial geometries: (i) a bottlenecked Josephson junction with external flux, (ii) a nanodisk ferromagnet deposited on top of a superconductor and (iii) superconducting islands in contact with a ferromagnet. In case (i), we demonstrate that one may control externally not only the geometrical array in which superconducting vortices arrange themselves, but also to cause coalescence and tune the number of vortices. In case (iii), we show that the supercurrent path can be tailored by incorporating magnetic elements in planar Josephson junctions which also lead to a strong modulation of the density of states. The finite element method presented herein paves the way for gaining insight in physical phenomena which have remained largely unexplored due to the complexity of solving the full quasiclassical equations in higher dimensions.Nonlinear differential equations (NLDEs) play a pivotal role in virtually all areas of physics. They are used to describe completely disparate phenomena ranging from the behavior of ocean waves to the elasticity of materials. Thus, techniques to solve such equations are of general interest as they provide a way to obtain insight in a number of different physical systems. NLDEs are known for being notoriously difficult to solve and, more often than not, a set of NLDEs describing a particular physical scenario has to be addressed as a distinct problem since general techniques to solve such equations are scarce.In quantum condensed matter physics, mesoscopic systems both in and out of equilibrium represent a very important arena where NLDEs are prevalent. A powerful tool used to describe such systems is the quasiclassical Keldysh theory, which has been reviewed in several works [1][2][3][4][5][6][7] . The theory is based on a Green function method which thus has a natural way of including disorder and other types of self-energies in the system. The quasiclassical Keldysh theory is capable of treating both ballistic systems and "dirty" systems. In the latter case, quasiparticles are elastically scattered within the mean free path l mfp causing the resulting motion to be diffusive. In essence, the quasiclassical theory is a perturbation expansion valid when all energy scales in the problem are much smaller than the Fermi energy E F . Conversely, all length scales in the system should be much larger than the Fermi wavelength. This situation is realized in a number of mesoscopic systems, including normal metals, superconductors and weakly polarized ferromagnets. Strongly polarized ferromagnets, where t...