From footing simulations to slope stability analyses, geomaterials are frequently described by non-associated elasto-plastic constitutive models. To ensure dissipativity and exact resolution of the constitutive model at each time step, an implicit formulation of the numerical method combined with the classical return mapping algorithm is often preferred. This work, in continuity with Acary et al. [1], proposes an alternative approach through the development of a monolithic solver capable of simultaneously solving the plastic constitutive equation, the potential hardening and the Coulomb friction with unilateral contact. This approach avoids in particular the challenging computation of the consistent tangent operator and enables the use of efficient and robust solvers that come from numerical optimization. In order to account for instances where, as it is the case with slope failures, the material is subjected to a loading phase followed by a flow phase exhibiting large deformations, we use the Material Point Method (MPM). The method employs both Eulerian and Lagrangian discretisation, which makes it well suited for this class of problems. The constitutive equations are written as a differential inclusion, inspired by the Implicit Standard Material (ISM) framework [2]. A semi-smooth Newton method and a non-smooth Gauss-Siedel method are used to solve the constitutive law within an implicit time scheme. The monolithic structure of the solver permits the implementation of non-conforming Dirichlet boundary conditions and interaction with rigid bodies. The relevance of the method is demonstrated with a footing problem and slope stability simulations where the existence of an analytical solution in limit analysis for perfect plasticity is established. The simulation demonstrates a high degree of correlation between the MPM and the analytical solution with regard to ultimate bearing capacity, as well as a good adequation between the MPM and the FEM stress and strain distributions.