Abstract Nonperturbative analytical and numerical methods are presented for the solution of the coupled nonlinear Maxwell–Schrödinger system of partial differential equations. These equations have been derived within the Hamiltonian or canonical formalism. The canonical approach to dynamics, which begins with the Maxwell and Schrödinger Lagrangians together with a Lorenz gauge fixing term, yields a set of first order Hamilton equations which form a well-posed initial value problem. That is, their solution is uniquely determined once the initial values for each of the dynamical variables are specified. The equations are closed since the Schrödinger wavefunction is chosen to be the source for the electromagnetic field and the electromagnetic field acts back upon the wavefunction, thus producing new fields. In practice, the Maxwell–Schrödinger Lagrangian is represented in a spatial basis of Gaussian functions. Application of the calculus of variations leads to a set of dynamical equations that, for that choice of basis, represent the coupled first order Maxwell–Schrödinger equations. In the limit of a complete basis these equations are exact and for any finite choice of basis they form an approximate system of dynamical equations that can be integrated in time and made systematically more accurate by enriching the basis. The dynamics of the basis-represented Maxwell–Schrödinger system is numerically investigated for a single spinless hydrogen atom interacting with the electromagnetic field in a cavity.