In many applications of practical interest, solutions of partial differential equation models arise as critical points of an underlying (energy) functional. If such solutions are saddle points, rather than being maxima or minima, then the theoretical framework is nonstandard, and the development of suitable numerical approximation procedures turns out to be highly challenging. In this paper, our aim is to present an iterative discretization methodology for the numerical solution of nonlinear variational problems with multiple (saddle point) solutions. In contrast to traditional numerical approximation schemes, which typically fail in such situations, the key idea of the current work is to employ a simultaneous interplay of a previously developed local minimax approach and adaptive Galerkin discretizations. We thereby derive an adaptive local minimax Galerkin (LMMG) method, which combines the search for saddle point solutions and their approximation in finite-dimensional spaces in a highly effective way. Under certain assumptions, we will prove that the generated sequence of approximate solutions converges to the solution set of the variational problem. This general framework will be applied to the specific context of finite element discretizations of (singularly perturbed) semilinear elliptic boundary value problems, and a series of numerical experiments will be presented.