To design and control magnetic confinement fusion reactors, one must compute the geometry of the confined plasma at equilibrium. In this talk, I present a new method for solving this free-boundary problem through integral equations and PDE-constrained optimization. I detail a formulation of the equilibrium equations that couples the Grad-Shafranov equation for the interior magnetic field with an integral equation for the exterior vacuum field. I introduce high order numerical methods for solving the associated integral equation, including a novel application of the Kapur-Rokhlin quadrature rule for singular integrals in axisymmetric magnetic confinement systems. Finally, I frame the coupled equations in the larger PDE-constrained optimization framework and discuss gradient descent methods for the optimization iteration.