A new perturbation and continuation method is presented for computing and analyzing stellarator equilibria. The method is formally derived from a series expansion about the equilibrium condition F ≡ J ×B−∇p = 0, and an efficient algorithm for computing solutions to 2nd and 3rd order perturbations is developed. The method has been implemented in the DESC stellarator equilibrium code, using automatic differentiation to compute the required derivatives. Examples are shown demonstrating its use for computing complicated equilibria, perturbing a tokamak into a stellarator, and performing parameter scans in pressure, rotational transform and boundary shape in a fraction of the time required for a full solution.