Many problems in engineering can be understood as controlling the bifurcation structure of a given device. For example, one may wish to delay the onset of instability, or bring forward a bifurcation to enable rapid switching between states. We propose a numerical technique for controlling the bifurcation diagram of a nonlinear partial differential equation by varying the shape of the domain. Specifically, we are able to delay or advance a given bifurcation point to a given parameter value, often to within machine precision. The algorithm consists of solving a shape optimization problem constrained by an augmented system of equations, the Moore-Spence system, that characterize the location of the bifurcation points. Numerical experiments on the Allen-Cahn, Navier-Stokes, and hyperelasticity equations demonstrate the effectiveness of this technique in a wide range of settings.