A fundamental task in many applications of modern scientific computing consists in the numerical solution of parametrized nonlinear systems \[F(\xb(\xib_j),\xib_j)=0,\quad j=1,\dots,N.\] When the number of parameter instances $N$ is large, solving these systems independently can be computationally prohibitive. While a rich literature exists on efficient algorithms for parametrized linear problems, the development of comparable methods for nonlinear systems remains far less mature. In this talk, we focus on a structured class of nonlinear systems and show that the associated snapshot matrix $X=[\xb(\xib_1)|\dots,|\xb(\xib_N)]$ exhibits a fast decay of singular values, indicating an underlying low-rank structure. To exploit this property, we reinterpret the collection of nonlinear systems as the first-order optimality conditions of an optimization problem posed over the set of real matrices. We thus compute a low-rank approximation to the snapshot matrix $X$ by minimizing a suitable functional over the manifold of fixed-rank matrices using Riemannian optimization methods, including Riemannian gradient and trust-region algorithms. To improve efficiency and robustness, suitable preconditioning strategies are employed, leading to robust convergence with respect to the size of the nonlinear systems and number of parameters. To handle the nonlinear contributions in a low-rank format, we leverage CUR decompositions based on DEIM approximations, together with rank adaptivity. Numerical experiments demonstrate the effectiveness of the proposed framework.