We present a continuous formulation of epidemic spreading on multilayer networks using a tensorial representation, extending the models of monoplex networks to this context. We derive analytical expressions for the epidemic threshold of the SIS and SIR dynamics, as well as upper and lower bounds for the disease prevalence in the steady state for the SIS scenario. Using the quasi-stationary state (QS) method we numerically show the existence of disease localization and the emergence of two or more phase transitions, which are characterized analytically and numerically through the inverse participation ratio. Furthermore, when mapping the critical dynamics to an eigenvalue problem, we observe a characteristic transition in the eigenvalue spectra of the supra-contact tensor as a function of the ratio of two spreading rates: If the rate at which the disease spreads within a layer is comparable to the spreading rate across layers, the individual spectra of each layer merge with the coupling between layers. The formalism introduced here provides a unifying mathematical approach to disease contagion in multiplex systems opening new possibilities for the study of spreading processes. Finally, our findings show the importance of considering the multilayer nature of many real systems, as this interdependency usually gives raise to new phenomenology.