The Einstein constraints equations are of fundamental interest in the the study of Einstein's theory of general relativity. This coupled nonlinear elliptic system must also be solved reliably and efficiently for gravitational wave simulation. The equations have been studied intensively for half a century; they are a particular example of a "critical exponent" problem, often arising in geometric analysis. In this lecture, we begin with an overview of the most useful mathematical formulation of the constraint equations, and then summarize the known existence, uniqueness, and multiplicity results through 2008. We then present a number of new existence and multiplicity results developed since 2008 that substantially extend the solution theory for the constraint equations. The techniques needed for developing new results are wide-ranging, and include fixed-point arguments, maximum principles, a priori estimates, bifurcation theory, and other techniques in nonlinear analysis and partial differential equations. We then shift gears a bit and consider Galerkin and Petrov-Galerkin type approximation methods for developing "good" numerical methods for solving this system. We examine how one proves rigorous error estimates for particular classes of numerical methods, including both classical finite element methods and newer methods from the finite element exterior calculus. This project is joint work with a number of collaborators over several years.