The contact of solids with rough surfaces plays a fundamental role in physical phenomena such as friction, wear, sealing, and thermal transfer. However, its simulation is a challenging problem due to surface asperities covering a wide range of length-scales. In addition, non-linear local processes, such as plasticity, are expected to occur even at the lightest loads. In this context, robust and efficient computational approaches are required. We therefore present a novel numerical method, based on integral equations, capable of handling the large discretization requirements of real rough surfaces as well as the non-linear plastic flow occurring below and at the contacting asperities. This method is based on a new derivation of the Mindlin fundamental solution in Fourier space, which leverages the computational efficiency of the fast Fourier transform. The use of this Mindlin solution allows a dramatic reduction of the memory imprint (as the Fourier coefficients are computed on-the-fly), a reduction of the discretization error, and the exploitation of the structure of the functions to speed up computation of the integral operators. We validate our method against an elasticplastic FEM Hertz normal contact simulation and showcase its ability to simulate contact of rough surfaces with plastic flow.