In this work, we further develop multigoal-oriented a posteriori error estimation with two objectives in mind. First, we formulate goal-oriented mesh adaptivity for multiple functionals of interest for nonlinear problems in which both the Partial Differential Equation (PDE) and the goal functionals may be nonlinear. Our method is based on a posteriori error estimates in which the adjoint problem is used and a partition-of-unity is employed for the error localization that allows us to formulate the error estimator in the weak form. We provide a careful derivation of the primal and adjoint parts of the error estimator. The second objective is concerned with balancing the nonlinear iteration error with the discretization error yielding adaptive stopping rules for Newton’s method. Our techniques are substantiated with several numerical examples including scalar PDEs and PDE systems, geometric singularities, and both nonlinear PDEs and nonlinear goal functionals. In these tests, up to six goal functionals are simultaneously controlled.