A mathematical model and numerical framework are presented for computing multi-physics multi-layer coating flow dynamics, with applications to the leveling of multi-layer paint films. The algorithm combines finite difference level set methods and high-order accurate sharp-interface implicit mesh discontinuous Galerkin methods to capture a complex set of multi-physics, incorporating Marangoni-driven multi-phase interfacial flow and the transport, mixing, and evaporation of multiple dissolved species. In particular, we develop several numerical methods for this multi-physics problem, including: high-order local discontinuous Galerkin methods for Poisson problems with Robin boundary conditions on implicitly-defined domains, to capture solvent evaporation; finite difference surface gradient methods, to robustly and accurately incorporate Marangoni stresses; and a coupled multi-physics time stepping approach, to incorporate all the different solvers at play including quasi-Newtonian fluid flow. The framework is applicable to an arbitrary number of layers and dissolved species; here, we apply it in a variety of settings, including multi-solvent evaporative paint dynamics, the flow and leveling of multi-layer automobile paint coatings in both 2D and 3D, and an examination of interfacial turbulence within a multi-layer matter cascade. Our results reproduce several phenomena observed in experiment, such as the formation of Marangoni plumes and Bénard cells. We also use the model to study the impact of long-wave deformational surface modes on immersed interfaces as well as the emergence of the final multi-layer film profile.