Nonlinear designs are notorious for their difficulty in getting rigorously treated. In this talk, the full electromagnetic interactions with nonlinear metasurfaces will be analytically described in several setups incorporating planar or cylindrical geometries, coupled or not, being excited normally or obliquely. The rigorous solutions are determined just by picking the roots of highly transcendental functions within a closed interval. Significant multistability, namely, the property of obtaining different outputs depending on the past values of the input which is the key effect behind memory components, has been exhibited by the devices. The reported results are expected to assist the modeling of nonlinear metasurfaces and open unexplored opportunities towards the efficient design of photonic memory elements.