Higher-order networks are gaining significant scientific attention due to their ability to encode the many-body interactions present in complex systems. However, higher-order networks have the limitation that they only capture many-body interactions of the same type. To address this limitation, we present a mathematical framework that determines the topology of higher-order multiplex networks and illustrates the interplay between their topology and dynamics. Specifically, we examine the diffusion of topological signals associated not only to the nodes but also to the links and to the higher-dimensional simplices of multiplex simplicial complexes. We leverage on the ubiquitous presence of the overlap of the simplices to couple the dynamics among multiplex layers, introducing a definition of multiplex Hodge Laplacians and Dirac operators. We show that the spectral properties of these operators determine higher-order diffusion on higher-order multiplex networks and encode their multiplex Betti numbers. Our numerical investigation of the spectral properties of synthetic and real (connectome, microbiome) multiplex simplicial complexes indicates that the coupling between the layers can either speed up or slow down the higher-order diffusion of topological signals. This mathematical framework is very general and can be applied to study generic higher-order systems with interactions of multiple types. In particular, these results might find applications in brain networks which are understood to be both multilayer and higher-order.