  We demonstrate for the first time that the nonlinear saturation of
neoclassical tearing modes (NTMs) can be found directly using a variational
principle based on Taylor relaxation, without needing to simulate the
intermediate, resistivity-dependent dynamics. As in previous investigations of
classical tearing mode saturation (Loizu et al. 2020; Loizu & Bonfiglio 2023),
we make use of SPEC (Hudson et al. 2012), an equilibrium solver based on the
variational principle of the Multi-Region relaxed MHD, featuring stepped
pressure profiles and arbitrary magnetic topology. We work in slab geometry and
employ a simple bootstrap current model $J_\textrm{bs} = C \nabla p$ to study
the bootstrap-driven tearing modes, scanning over the asymptotic matching
parameter $\Delta'$ and the bootstrap current strength. Saturated island widths
produced by SPEC agree well with the predictions of an initial value resistive
MHD code (Huang & Bhattacharjee 2016) while being orders of magnitude faster to
calculate. Additionally, we observe good agreement with a simple analytical
Modified Rutherford Equation, without requiring any fitting coefficients. The
match is obtained for both linearly unstable classical tearing modes in the
presence of bootstrap current, and neoclassical tearing modes, which are
linearly stable but nonlinear-unstable due to the effects of the bootstrap
current
