degradation functions in which the derivative $g^\prime \left( \phi \right)$ is nonlinear, the subsystem represented by also becomes nonlinear, resulting in the need to perform nested iterations. In our experience, a naive implementation of the alternate minimization algorithm wherein the linearized phase-field equation is solved only once before going back to linear momentum is generally unstable and leads to incorrect results. That is, in the subsystem $$\left\{ \phi \right\}_i^{m+1} = \left\{ \phi \right\}_i^m - \left[ \bm{K}^{\phi\phi} \left( \phi_i^m \right) \right]^{-1} \left\{ r^{\phi} \right\}_i^m$$ corresponding to the $\left( m+1 \right)$th iteration within the $i$th time step, use of the exact Jacobian given by $$\bm{K}^{\phi\phi}_{IJ} \left( \phi_i^m \right) = \int_\Omega \left[ \mathcal{G}_c\ell\, \bm{B}_I^{\phi T} \bm{B}_J^\phi + \left( \frac{\mathcal{G}_c}{\ell} + \mathcal{H} g^{\prime\prime} \left( \phi_i^m \right) \right) N_I N_J \right] {\,\mathrm{d}}\Omega$$ often produces an incorrect evolution of the phase-field and eventual blow-up. However, we have found that such an approach can be made to work by replacing $g^{\prime\prime} \left( \phi_i^m \right)$ with an approximate expression derived from low order terms in a Taylor expansion. Assuming that the degradation function is smooth, we can obtain $g^\prime \left( \phi_1 \right)$ as an infinite sum of terms involving higher order derivatives of $g$ evaluated at $\phi_2 \in 