  We present a boundary integral equation solver for computing Taylor relaxed
states in non-axisymmetric solid and shell-like toroidal geometries. The
computation of Taylor states in these geometries is a key element for the
calculation of stepped pressure stellarator equilibria. The integral
representation of the magnetic field in this work is based on the generalized
Debye source formulation, and results in a well-conditioned second-kind
boundary integral equation. The integral equation solver is based on a spectral
discretization of the geometry and unknowns, and the computation of the
associated weakly-singular integrals is performed with high-order quadrature
based on a partition of unity. The resulting scheme for applying the integral
operator is then coupled with an iterative solver and suitable preconditioners.
Several numerical examples are provided to demonstrate the accuracy and
efficiency of our method, and a direct comparison with the leading code in the
field is reported.
