diff --git a/.azure/templates/install.yml b/.azure/templates/install.yml
index 7b85780efe..6c0283bea6 100644
--- a/.azure/templates/install.yml
+++ b/.azure/templates/install.yml
@@ -53,9 +53,11 @@ steps:
     displayName: Install Anaconda packages
 
   - bash: |
+      source activate myEnvironment
       if [[ ${{ parameters.python_version }} != "2.7" ]]; then
-        source activate myEnvironment
         pip install scikit-fmm
         pip install packaging
       fi
+      pip install git+https://github.com/guyer/steppyngstounes.git
     displayName: Install pip packages
+
diff --git a/documentation/USAGE.rst b/documentation/USAGE.rst
index 8ce1d71ca6..49d7f24220 100644
--- a/documentation/USAGE.rst
+++ b/documentation/USAGE.rst
@@ -1050,6 +1050,19 @@ Nothing different needs to be done when
 
    %    \subsection{Internal boundary conditions}
 
+-----------------
+Adaptive Stepping
+-----------------
+
+Step size can be controlled with the :term:`steppyngstounes` package.
+Demonstrations of its use are found in :mod:`examples.phase.binary` and
+:mod:`examples.phase.binaryCoupled.`
+
+.. note::
+
+   The old :mod:`fipy.steppers` classes are now deprecated.  They were
+   undocumented and did not work very well.
+
 .. _RunningUnderPython2:
 
 ----------------------
diff --git a/documentation/glossary.rst b/documentation/glossary.rst
index 346d64e966..6c857d3a88 100644
--- a/documentation/glossary.rst
+++ b/documentation/glossary.rst
@@ -195,6 +195,11 @@ Glossary
       See
       http://sphinx.pocoo.org/.
 
+   steppyngstounes
+      This package provides iterators that simplify both deterministic and
+      adaptive time (or other independent variable) stepping.
+      See https://github.com/guyer/steppyngstounes.
+
    TravisCI
       A cloud-based :term:`Continuous Integration` tool.
       See https://travis-ci.org.
diff --git a/examples/phase/binary.py b/examples/phase/binary.py
index f3938969a2..2b8f0456c9 100755
--- a/examples/phase/binary.py
+++ b/examples/phase/binary.py
@@ -461,12 +461,70 @@ def deltaChemPot(phase, C, T):
 >>> sharp.setValue(Cs, where=x < L * fraction)
 >>> sharp.setValue(Cl, where=x >= L * fraction)
 
+>>> elapsed = Variable(value=0.) # s
+
 .. index::
    module: fipy.viewers
 
 >>> if __name__ == '__main__':
-...     viewer = Viewer(vars=(phase, C, sharp),
-...                     datamin=0., datamax=1.)
+...     try:
+...         from fipy import Matplotlib1DViewer
+...         from scipy.optimize import leastsq # doctest: +SCIPY
+...         from matplotlib import pyplot as plt
+...
+...         class PhaseViewer(Matplotlib1DViewer):
+...             def __init__(self, phase, C, sharp, elapsed, title=None,
+...                          tmin=None, tmax=None, **kwlimits):
+...                 self.phase = phase
+...                 self.x = self.phase.mesh.cellCenters[0]
+...                 self.elapsed = elapsed
+...
+...                 fig = plt.figure(constrained_layout=True, figsize=[4, 6])
+...                 gs = fig.add_gridspec(ncols=1, nrows=3)
+...                 self.viewax = fig.add_subplot(gs[1:, 0])
+...                 self.viewax.set_xlabel(r"$x / \mathrm{cm}$")
+...                 self.posax = fig.add_subplot(gs[0, 0], sharex=self.viewax)
+...                 self.posax.set_ylabel(r"$t / \mathrm{s}$")
+...                 self.posax.label_outer()
+...                 self.posax.set_ylim(tmin, tmax)
+...
+...                 self.times = []
+...                 self.positions = []
+...                 self.line, = self.posax.semilogy([0.], [0.],
+...                                                  color="blue")
+...                 self.bug, = self.posax.semilogy([], [],
+...                                                 linestyle="", marker="o", color="red")
+...
+...                 super(PhaseViewer, self).__init__(vars=(phase, C, sharp),
+...                                                   axes=self.viewax,
+...                                                   **kwlimits)
+...
+...                 self.timer = self.viewax.text(0.00025, 0.2, "t = 0")
+...
+...             @staticmethod
+...             def tanhResiduals(p, y, x):
+...                 x0, d = p
+...                 return y - 0.5 * (1 - numerix.tanh((x - x0) / (2*d)))
+...
+...             def _plot(self):
+...                 (x0_fit, d_fit), msg = leastsq(self.tanhResiduals, [L/2., deltaA],
+...                                                args=(self.phase.globalValue,
+...                                                      self.x.globalValue))
+...                 self.positions.append(x0_fit)
+...                 self.times.append(float(self.elapsed))
+...                 self.line.set_data(self.positions, self.times)
+...                 self.bug.set_data([self.positions[-1]], [self.times[-1]])
+...
+...                 self.timer.set_text(r"$t = {:f}\\,\\mathrm{{s}}$".format(float(self.elapsed)))
+...
+...                 super(PhaseViewer, self)._plot()
+...
+...         viewer = PhaseViewer(phase=phase, C=C, sharp=sharp,
+...                              elapsed=elapsed, tmin=1e-5, tmax=300 * 3600,
+...                              datamin=0., datamax=1.)
+...     except ImportError:
+...         viewer = Viewer(vars=(phase, C, sharp),
+...                         datamin=0., datamax=1.)
 ...     viewer.plot()
 
 Because the phase field interface will not move, and because we've seen in
@@ -496,7 +554,7 @@ def deltaChemPot(phase, C, T):
 >>> C.updateOld()
 >>> phaseRes = 1e+10
 >>> diffRes = 1e+10
->>> while phaseRes > 1e-3 or diffRes > 1e-3:
+>>> while phaseRes > 1e-3 or diffRes > 1e-3 or abs(Cavg.value - 0.5) > 5e-7:
 ...     phaseRes = phaseEq.sweep(var=phase, dt=dt)
 ...     diffRes = diffusionEq.sweep(var=C, dt=dt, solver=solver)
 >>> from fipy import input
@@ -505,7 +563,7 @@ def deltaChemPot(phase, C, T):
 ...     input("Stationary phase field. Press <return> to proceed...")
 
 .. image:: binary/stationary.*
-   :width: 90%
+   :width: 50%
    :align: center
    :alt: phase and composition fields in equilibrium, compared with phase diagram concentrations
 
@@ -553,14 +611,14 @@ def deltaChemPot(phase, C, T):
    \\vec{u}_\\phi &= \\frac{D_\\phi}{C} \\nabla \\phi
    \\\\
    &\\approx
-   \\frac{Dl \\frac{1}{2} V_m}{R T}
+   \\frac{D_l \\frac{1}{2} V_m}{R T}
    \\left[
        \\frac{L_B\\left(T - T_M^B\\right)}{T_M^B}
        - \\frac{L_A\\left(T - T_M^A\\right)}{T_M^A}
    \\right] \\frac{1}{\\Delta x}
    \\\\
    &\\approx
-   \\frac{Dl \\frac{1}{2} V_m}{R T}
+   \\frac{D_l \\frac{1}{2} V_m}{R T}
    \\left(L_B + L_A\\right) \\frac{T_M^A - T_M^B}{T_M^A + T_M^B}
    \\frac{1}{\\Delta x}
    \\\\
@@ -569,22 +627,18 @@ def deltaChemPot(phase, C, T):
 To get a :math:`\\text{CFL} = \\vec{u}_\\phi \\Delta t / \\Delta x < 1`, we need a
 time step of about :math:`\\unit{10^{-5}}{\\second}`.
 
->>> dt = 1.e-5
-
->>> if __name__ == '__main__':
-...     timesteps = 100
-... else:
-...     timesteps = 10
+>>> dt0 = 1.e-5
 
 >>> from builtins import range
->>> for i in range(timesteps):
+>>> for i in range(8):
 ...     phase.updateOld()
 ...     C.updateOld()
 ...     phaseRes = 1e+10
 ...     diffRes = 1e+10
-...     while phaseRes > 1e-3 or diffRes > 1e-3:
-...         phaseRes = phaseEq.sweep(var=phase, dt=dt)
-...         diffRes = diffusionEq.sweep(var=C, dt=dt, solver=solver)
+...     while phaseRes > 1e-3 or diffRes > 1e-3 or abs(Cavg.value - 0.5) > 1e-6:
+...         phaseRes = phaseEq.sweep(var=phase, dt=dt0)
+...         diffRes = diffusionEq.sweep(var=C, dt=dt0, solver=solver)
+...     elapsed.value = (i + 1) * dt0
 ...     if __name__ == '__main__':
 ...         viewer.plot()
 
@@ -592,17 +646,111 @@ def deltaChemPot(phase, C, T):
 >>> if __name__ == '__main__':
 ...     input("Moving phase field. Press <return> to proceed...")
 
-.. image:: binary/moving.*
-   :width: 90%
-   :align: center
-   :alt: phase and composition fields in during solidification, compared with final phase diagram concentrations
-
-We see that the composition on either side of the interface approach the
+We see that the composition on either side of the interface approaches the
 sharp-interface solidus and liquidus, but it will take a great many more
 timesteps to reach equilibrium. If we waited sufficiently long, we
 could again verify the final concentrations and phase fraction against the
 expected values.
 
+We can estimate the time to equilibration by examining the time for the
+diffusion field to become uniform.  In the liquid, this will take
+:math:`\\mathcal{O}((\\unit{10}{\\micro\\meter})^2 / D_l) =
+\\unit{0.1}{\\second}` and in the solid
+:math:`\\mathcal{O}((\\unit{10}{\\micro\\meter})^2 / D_s) =
+\\unit{1000}{\\second}`.
+
+Not wanting to take a hundred-million steps, we employ adaptive time
+stepping, using the :term:`steppyingstounes` package.  This package takes
+care of many of the messy details of stepping, like overshoot, underflow,
+and step size adaptation, while keeping the structure of our solve loop
+largely intact.
+
+>>> from steppyngstounes import SequenceStepper, PIDStepper
+>>> from itertools import count
+
+Assuming the process is dominated by diffusion, we can take steps that
+increase geometrically.  Since we're unsure if diffusion is the only
+process controlling dynamics, we take each increasing step with an adaptive
+stepper that uses a `PID controller`_ to keep the equation residuals and
+mass conservation within acceptable limits.  The total number of solves is
+not strongly sensitive to the number of sweeps, but two sweeps seems to be
+both sufficient and efficient.
+
+We'll only advance the step if it's successful, so we need to update the
+old values before we get started.
+
+>>> phase.updateOld()
+>>> C.updateOld()
+
+>>> if __name__ == '__main__':
+...     totaltime = 300 * 3600 # 300 h
+... else:
+...     totaltime = 32e-5 # 320 us
+
+>>> dt = dt0
+
+>>> for checkpoint in SequenceStepper(start=float(elapsed), stop=totaltime,
+...                                   sizes=(dt0 * 2**(n/2) for n in count(7))):
+...     for step in PIDStepper(start=checkpoint.begin,
+...                            stop=checkpoint.end,
+...                            size=dt):
+...         for sweep in range(2):
+...             phaseRes = phaseEq.sweep(var=phase, dt=step.size)
+...             diffRes = diffusionEq.sweep(var=C, dt=step.size, solver=solver)
+...         err = max(phaseRes / 1e-3,
+...                   diffRes / 1e-3,
+...                   abs(Cavg.value - 0.5) / 1e-6)
+...         if step.succeeded(error=err):
+...             phase.updateOld()
+...             C.updateOld()
+...             elapsed.value = step.end
+...         else:
+...             phase.value = phase.old
+...             C.value = C.old
+...     # the last step might have been smaller than possible,
+...     # if it was near the end of the checkpoint range
+...     dt = step.want
+...     _ = checkpoint.succeeded()
+...     if __name__ == '__main__':
+...         viewer.plot()
+
+>>> from fipy import input
+>>> if __name__ == '__main__':
+...     input("Re-equilbrated phase field. Press <return> to proceed...")
+
+.. image:: binary/binary-0.000899.*
+   :width: 30%
+   :alt: phase and composition fields at t=0.000899, compared with final phase diagram concentrations
+
+.. image:: binary/binary-8.949963.*
+   :width: 30%
+   :alt: phase and composition fields at t=8.949963, compared with final phase diagram concentrations
+
+.. image:: binary/binary-1080000.000000.*
+   :width: 30%
+   :alt: phase and composition fields at t=1080000, compared with final phase diagram concentrations
+
+The interface moves :math:`\\approx \\unit{3.4}{\\micro\\meter}` in
+:math:`\\unit{80}{\\milli\\second}`, driven by diffusion in the liquid
+phase (compare the estimate above of :math:`\\unit{0.1}{\\second}`).
+For the next
+:math:`\\unit{20}{\\second}`, the interface stalls while the solute step
+trapped in the solid phase diffuses outward
+(:math:`(\\unit{3.4}{\\micro\\meter})^2 / D_s =
+\mathcal{O}(\\unit{100}{\\second})`).  Once the solute gradient in the
+solid reaches the new position of the interface, the solidification front
+begins to move, driven by diffusion in the solid.  When the solute in the
+solid becomes uniform, the interface stalls again after :math:`\\approx
+\\unit{4000}{\\second}`, having moved another
+:math:`\\unit{3.2}{\\micro\\meter}` (recall the estimate of
+:math:`\\unit{1000}{\\second}` for equilibration in the solid).  After this
+point, there is essentially no further motion of the interface and barely
+perceptible changes in the concentration field.  The fact that the
+interface does not reach the predicted phase fraction is due to the fact
+that this phase field model accounts for adsorption at the interface,
+resulting in the bulk phases not having exactly the concentrations assumed
+in the sharp interface treatment.
+
 .. rubric:: Footnotes
 
 .. [#phi] We will find that we need to "sweep" this non-linear problem
@@ -615,6 +763,8 @@ def deltaChemPot(phase, C, T):
    as a :class:`~fipy.variables.variable.Variable`
 
 .. _CFL limit: http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition
+
+.. _PID controller: https://en.wikipedia.org/wiki/PID_controller
 """
 
 from __future__ import unicode_literals
diff --git a/examples/phase/binaryCoupled.py b/examples/phase/binaryCoupled.py
index 6324ce3fa6..895ecbc207 100755
--- a/examples/phase/binaryCoupled.py
+++ b/examples/phase/binaryCoupled.py
@@ -479,12 +479,70 @@ def deltaChemPot(phase, C, T):
 >>> sharp.setValue(Cs, where=x < L * fraction)
 >>> sharp.setValue(Cl, where=x >= L * fraction)
 
+>>> elapsed = Variable(value=0.) # s
+
 .. index::
    module: fipy.viewers
 
 >>> if __name__ == '__main__':
-...     viewer = Viewer(vars=(phase, C, sharp),
-...                     datamin=0., datamax=1.)
+...     try:
+...         from fipy import Matplotlib1DViewer
+...         from scipy.optimize import leastsq # doctest: +SCIPY
+...         from matplotlib import pyplot as plt
+...
+...         class PhaseViewer(Matplotlib1DViewer):
+...             def __init__(self, phase, C, sharp, elapsed, title=None,
+...                          tmin=None, tmax=None, **kwlimits):
+...                 self.phase = phase
+...                 self.x = self.phase.mesh.cellCenters[0]
+...                 self.elapsed = elapsed
+...
+...                 fig = plt.figure(constrained_layout=True, figsize=[4, 6])
+...                 gs = fig.add_gridspec(ncols=1, nrows=3)
+...                 self.viewax = fig.add_subplot(gs[1:, 0])
+...                 self.viewax.set_xlabel(r"$x / \mathrm{cm}$")
+...                 self.posax = fig.add_subplot(gs[0, 0], sharex=self.viewax)
+...                 self.posax.set_ylabel(r"$t / \mathrm{s}$")
+...                 self.posax.label_outer()
+...                 self.posax.set_ylim(tmin, tmax)
+...
+...                 self.times = []
+...                 self.positions = []
+...                 self.line, = self.posax.semilogy([0.], [0.],
+...                                                  color="blue")
+...                 self.bug, = self.posax.semilogy([], [],
+...                                                 linestyle="", marker="o", color="red")
+...
+...                 super(PhaseViewer, self).__init__(vars=(phase, C, sharp),
+...                                                   axes=self.viewax,
+...                                                   **kwlimits)
+...
+...                 self.timer = self.viewax.text(0.00025, 0.2, "t = 0")
+...
+...             @staticmethod
+...             def tanhResiduals(p, y, x):
+...                 x0, d = p
+...                 return y - 0.5 * (1 - numerix.tanh((x - x0) / (2*d)))
+...
+...             def _plot(self):
+...                 (x0_fit, d_fit), msg = leastsq(self.tanhResiduals, [L/2., deltaA],
+...                                                args=(self.phase.globalValue,
+...                                                      self.x.globalValue))
+...                 self.positions.append(x0_fit)
+...                 self.times.append(float(self.elapsed))
+...                 self.line.set_data(self.positions, self.times)
+...                 self.bug.set_data([self.positions[-1]], [self.times[-1]])
+...
+...                 self.timer.set_text(r"$t = {:f}\\,\\mathrm{{s}}$".format(float(self.elapsed)))
+...
+...                 super(PhaseViewer, self)._plot()
+...
+...         viewer = PhaseViewer(phase=phase, C=C, sharp=sharp,
+...                              elapsed=elapsed, tmin=1e-5, tmax=300 * 3600,
+...                              datamin=0., datamax=1.)
+...     except ImportError:
+...         viewer = Viewer(vars=(phase, C, sharp),
+...                         datamin=0., datamax=1.)
 ...     viewer.plot()
 
 Because the phase field interface will not move, and because we've seen in
@@ -516,7 +574,7 @@ def deltaChemPot(phase, C, T):
 >>> initialRes = None
 >>> sweep = 0
 
->>> while res > 1e-4 and sweep < 20:
+>>> while res > 1e-8 and sweep < 100:
 ...     res = eq.sweep(dt=dt, solver=solver)
 ...     if initialRes is None:
 ...         initialRes = res
@@ -528,9 +586,10 @@ def deltaChemPot(phase, C, T):
 ...     viewer.plot()
 ...     input("Stationary phase field. Press <return> to proceed...")
 
-.. image:: binary/stationary.*
-   :width: 90%
+.. image:: binary/coupled-stationary.*
+   :width: 50%
    :align: center
+   :alt: phase and composition fields in equilibrium, compared with phase diagram concentrations
 
 We verify that the bulk phases have shifted to the predicted solidus and
 liquidus compositions
@@ -578,14 +637,14 @@ def deltaChemPot(phase, C, T):
    \\vec{u}_\\phi &= \\frac{D_\\phi}{C} \\nabla \\phi
    \\\\
    &\\approx
-   \\frac{Dl \\frac{1}{2} V_m}{R T}
+   \\frac{D_l \\frac{1}{2} V_m}{R T}
    \\left[
        \\frac{L_B\\left(T - T_M^B\\right)}{T_M^B}
        - \\frac{L_A\\left(T - T_M^A\\right)}{T_M^A}
    \\right] \\frac{1}{\\Delta x}
    \\\\
    &\\approx
-   \\frac{Dl \\frac{1}{2} V_m}{R T}
+   \\frac{D_l \\frac{1}{2} V_m}{R T}
    \\left(L_B + L_A\\right) \\frac{T_M^A - T_M^B}{T_M^A + T_M^B}
    \\frac{1}{\\Delta x}
    \\\\
@@ -594,22 +653,18 @@ def deltaChemPot(phase, C, T):
 To get a :math:`\\text{CFL} = \\vec{u}_\\phi \\Delta t / \\Delta x < 1`, we need a
 time step of about :math:`\\unit{10^{-5}}{\\second}`.
 
->>> dt = 1.e-5
-
->>> if __name__ == '__main__':
-...     timesteps = 100
-... else:
-...     timesteps = 10
+>>> dt0 = 1.e-5
 
 >>> from builtins import range
->>> for i in range(timesteps):
+>>> for i in range(8):
 ...     phase.updateOld()
 ...     C.updateOld()
 ...     res = 1e+10
 ...     sweep = 0
-...     while res > 1e-3 and sweep < 20:
-...         res = eq.sweep(dt=dt, solver=solver)
+...     while (res > 1e-3 or abs(Cavg.value - 0.5) > 1e-8) and sweep < 20:
+...         res = eq.sweep(dt=dt0, solver=solver)
 ...         sweep += 1
+...     elapsed.value = (i + 1) * dt0
 ...     if __name__ == '__main__':
 ...         viewer.plot()
 
@@ -617,16 +672,112 @@ def deltaChemPot(phase, C, T):
 >>> if __name__ == '__main__':
 ...     input("Moving phase field. Press <return> to proceed...")
 
-.. image:: binary/moving.*
-   :width: 90%
-   :align: center
-
-We see that the composition on either side of the interface approach the
+We see that the composition on either side of the interface approaches the
 sharp-interface solidus and liquidus, but it will take a great many more
 timesteps to reach equilibrium. If we waited sufficiently long, we
 could again verify the final concentrations and phase fraction against the
 expected values.
 
+We can estimate the time to equilibration by examining the time for the
+diffusion field to become uniform.  In the liquid, this will take
+:math:`\\mathcal{O}((\\unit{10}{\\micro\\meter})^2 / D_l) =
+\\unit{0.1}{\\second}` and in the solid
+:math:`\\mathcal{O}((\\unit{10}{\\micro\\meter})^2 / D_s) =
+\\unit{1000}{\\second}`.
+
+Not wanting to take a hundred-million steps, we employ adaptive time
+stepping, using the :term:`steppyingstounes` package.  This package takes
+care of many of the messy details of stepping, like overshoot, underflow,
+and step size adaptation, while keeping the structure of our solve loop
+largely intact.
+
+>>> from steppyngstounes import SequenceStepper, PIDStepper
+>>> from itertools import count
+
+Assuming the process is dominated by diffusion, we can take steps that
+increase geometrically.  Since we're unsure if diffusion is the only
+process controlling dynamics, we take each increasing step with an adaptive
+stepper that uses a `PID controller`_ to keep the equation residuals and
+mass conservation within acceptable limits.  The total number of solves is
+not strongly sensitive to the number of sweeps, but two sweeps seems to be
+both sufficient and efficient.
+
+We'll only advance the step if it's successful, so we need to update the
+old values before we get started.
+
+>>> phase.updateOld()
+>>> C.updateOld()
+
+>>> if __name__ == '__main__':
+...     totaltime = 300 * 3600 # 300 h
+... else:
+...     totaltime = 32e-5 # 320 us
+
+>>> dt = dt0
+
+>>> for checkpoint in SequenceStepper(start=float(elapsed), stop=totaltime,
+...                                   sizes=(dt0 * 2**(n/2) for n in count(7))):
+...     for step in PIDStepper(start=checkpoint.begin,
+...                            stop=checkpoint.end,
+...                            size=dt):
+...         for sweep in range(2):
+...             res = eq.sweep(dt=step.size, solver=solver)
+...         err = max(res / 1e-3,
+...                   abs(Cavg.value - 0.5) / 1e-6)
+...         if step.succeeded(error=err):
+...             phase.updateOld()
+...             C.updateOld()
+...             elapsed.value = step.end
+...         else:
+...             phase.value = phase.old
+...             C.value = C.old
+...     # the last step might have been smaller than possible,
+...     # if it was near the end of the checkpoint range
+...     dt = step.want
+...     _ = checkpoint.succeeded()
+...     if __name__ == '__main__':
+...         viewer.plot()
+
+>>> from fipy import input
+>>> if __name__ == '__main__':
+...     input("Re-equilbrated phase field. Press <return> to proceed...")
+
+.. image:: binary/binaryCoupled-0.000899.*
+   :width: 30%
+   :alt: phase and composition fields at t=0.000899, compared with final phase diagram concentrations
+
+.. image:: binary/binaryCoupled-8.949963.*
+   :width: 30%
+   :alt: phase and composition fields at t=8.949963, compared with final phase diagram concentrations
+
+.. image:: binary/binaryCoupled-1080000.000000.*
+   :width: 30%
+   :alt: phase and composition fields at t=1080000, compared with final phase diagram concentrations
+
+The interface moves :math:`\\approx \\unit{2.8}{\\micro\\meter}` in
+:math:`\\unit{70}{\\milli\\second}`, driven by diffusion in the liquid
+phase (compare the estimate above of :math:`\\unit{0.1}{\\second}`).
+For the next
+:math:`\\unit{12}{\\second}`, the interface stalls while the solute step
+trapped in the solid phase diffuses outward
+(:math:`(\\unit{2.8}{\\micro\\meter})^2 / D_s =
+\mathcal{O}(\\unit{80}{\\second})`).  Once the solute gradient in the
+solid reaches the new position of the interface, the solidification front
+begins to move, driven by diffusion in the solid.  When the solute in the
+solid becomes uniform, the interface stalls again after :math:`\\approx
+\\unit{4000}{\\second}`, having moved another
+:math:`\\unit{2.9}{\\micro\\meter}` (recall the estimate of
+:math:`\\unit{1000}{\\second}` for equilibration in the solid).  After this
+point, there is essentially no further motion of the interface and barely
+perceptible changes in the concentration field.
+
+.. note::
+
+   This evolution is qualitatively consistent with that seen in
+   :mod:`examples.phase.binary`, but the interface does not move as far and
+   the bulk concentrations are further from the phase diagram values.  The
+   computation also takes substantially longer than the uncoupled variant.
+
 .. rubric:: Footnotes
 
 .. [#phi] We will find that we need to "sweep" this non-linear problem
@@ -639,6 +790,8 @@ def deltaChemPot(phase, C, T):
    as a :class:`~fipy.variables.variable.Variable`
 
 .. _CFL limit: http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition
+
+.. _PID controller: https://en.wikipedia.org/wiki/PID_controller
 """
 from __future__ import unicode_literals
 
diff --git a/examples/phase/generated/binary/binary-0.000899.pdf b/examples/phase/generated/binary/binary-0.000899.pdf
new file mode 100644
index 0000000000..63a8fa3940
Binary files /dev/null and b/examples/phase/generated/binary/binary-0.000899.pdf differ
diff --git a/examples/phase/generated/binary/binary-0.000899.png b/examples/phase/generated/binary/binary-0.000899.png
new file mode 100644
index 0000000000..6d2c5722b4
Binary files /dev/null and b/examples/phase/generated/binary/binary-0.000899.png differ
diff --git a/examples/phase/generated/binary/binary-1080000.000000.pdf b/examples/phase/generated/binary/binary-1080000.000000.pdf
new file mode 100644
index 0000000000..336b628264
Binary files /dev/null and b/examples/phase/generated/binary/binary-1080000.000000.pdf differ
diff --git a/examples/phase/generated/binary/binary-1080000.000000.png b/examples/phase/generated/binary/binary-1080000.000000.png
new file mode 100644
index 0000000000..36587e81c3
Binary files /dev/null and b/examples/phase/generated/binary/binary-1080000.000000.png differ
diff --git a/examples/phase/generated/binary/binary-8.949963.pdf b/examples/phase/generated/binary/binary-8.949963.pdf
new file mode 100644
index 0000000000..a2c7ffc216
Binary files /dev/null and b/examples/phase/generated/binary/binary-8.949963.pdf differ
diff --git a/examples/phase/generated/binary/binary-8.949963.png b/examples/phase/generated/binary/binary-8.949963.png
new file mode 100644
index 0000000000..0f4fa8b1c4
Binary files /dev/null and b/examples/phase/generated/binary/binary-8.949963.png differ
diff --git a/examples/phase/generated/binary/binaryCoupled-0.000899.pdf b/examples/phase/generated/binary/binaryCoupled-0.000899.pdf
new file mode 100644
index 0000000000..6bd3fd49d0
Binary files /dev/null and b/examples/phase/generated/binary/binaryCoupled-0.000899.pdf differ
diff --git a/examples/phase/generated/binary/binaryCoupled-0.000899.png b/examples/phase/generated/binary/binaryCoupled-0.000899.png
new file mode 100644
index 0000000000..c74fd2ab53
Binary files /dev/null and b/examples/phase/generated/binary/binaryCoupled-0.000899.png differ
diff --git a/examples/phase/generated/binary/binaryCoupled-1080000.000000.pdf b/examples/phase/generated/binary/binaryCoupled-1080000.000000.pdf
new file mode 100644
index 0000000000..e15c0e7fe3
Binary files /dev/null and b/examples/phase/generated/binary/binaryCoupled-1080000.000000.pdf differ
diff --git a/examples/phase/generated/binary/binaryCoupled-1080000.000000.png b/examples/phase/generated/binary/binaryCoupled-1080000.000000.png
new file mode 100644
index 0000000000..5994c798b1
Binary files /dev/null and b/examples/phase/generated/binary/binaryCoupled-1080000.000000.png differ
diff --git a/examples/phase/generated/binary/binaryCoupled-8.949963.pdf b/examples/phase/generated/binary/binaryCoupled-8.949963.pdf
new file mode 100644
index 0000000000..6835e184c0
Binary files /dev/null and b/examples/phase/generated/binary/binaryCoupled-8.949963.pdf differ
diff --git a/examples/phase/generated/binary/binaryCoupled-8.949963.png b/examples/phase/generated/binary/binaryCoupled-8.949963.png
new file mode 100644
index 0000000000..80bb0f69a7
Binary files /dev/null and b/examples/phase/generated/binary/binaryCoupled-8.949963.png differ
diff --git a/examples/phase/generated/binary/coupled-stationary.pdf b/examples/phase/generated/binary/coupled-stationary.pdf
new file mode 100644
index 0000000000..02edc6aab6
Binary files /dev/null and b/examples/phase/generated/binary/coupled-stationary.pdf differ
diff --git a/examples/phase/generated/binary/coupled-stationary.png b/examples/phase/generated/binary/coupled-stationary.png
new file mode 100644
index 0000000000..6a194f61c9
Binary files /dev/null and b/examples/phase/generated/binary/coupled-stationary.png differ
diff --git a/examples/phase/generated/binary/moving.pdf b/examples/phase/generated/binary/moving.pdf
deleted file mode 100644
index 422f7f6232..0000000000
Binary files a/examples/phase/generated/binary/moving.pdf and /dev/null differ
diff --git a/examples/phase/generated/binary/moving.png b/examples/phase/generated/binary/moving.png
deleted file mode 100644
index c1f11a33b7..0000000000
Binary files a/examples/phase/generated/binary/moving.png and /dev/null differ
diff --git a/examples/phase/generated/binary/stationary.pdf b/examples/phase/generated/binary/stationary.pdf
index 91cce042d8..b93bbd800a 100644
Binary files a/examples/phase/generated/binary/stationary.pdf and b/examples/phase/generated/binary/stationary.pdf differ
diff --git a/examples/phase/generated/binary/stationary.png b/examples/phase/generated/binary/stationary.png
index 099453ac24..9c8e29d0bc 100644
Binary files a/examples/phase/generated/binary/stationary.png and b/examples/phase/generated/binary/stationary.png differ
diff --git a/examples/phase/index.rst b/examples/phase/index.rst
index e6bae7acae..ea9d83ada4 100644
--- a/examples/phase/index.rst
+++ b/examples/phase/index.rst
@@ -7,6 +7,7 @@ Phase Field Examples
    :template: example.rst
 
    examples.phase.simple
+   examples.phase.binary
    examples.phase.binaryCoupled
    examples.phase.quaternary
    examples.phase.anisotropy
diff --git a/fipy/__init__.py b/fipy/__init__.py
index dce6a9797c..a65825e4d5 100644
--- a/fipy/__init__.py
+++ b/fipy/__init__.py
@@ -66,7 +66,6 @@ def excepthook(*args):
 from fipy.boundaryConditions import *
 from fipy.meshes import *
 from fipy.solvers import *
-from fipy.steppers import *
 from fipy.terms import *
 from fipy.tools import *
 from fipy.variables import *
@@ -76,7 +75,6 @@ def excepthook(*args):
 __all__.extend(boundaryConditions.__all__)
 __all__.extend(meshes.__all__)
 __all__.extend(solvers.__all__)
-__all__.extend(steppers.__all__)
 __all__.extend(terms.__all__)
 __all__.extend(tools.__all__)
 __all__.extend(variables.__all__)
diff --git a/fipy/steppers/__init__.py b/fipy/steppers/__init__.py
index 360029815f..701f59dd42 100644
--- a/fipy/steppers/__init__.py
+++ b/fipy/steppers/__init__.py
@@ -5,11 +5,14 @@
 from fipy.steppers.stepper import Stepper
 from fipy.steppers.pseudoRKQSStepper import PseudoRKQSStepper
 from fipy.steppers.pidStepper import PIDStepper
+from fipy.tools.decorators import deprecate
 
 __all__ = ["L1error", "L2error", "LINFerror", "sweepMonotonic"]
 from future.utils import text_to_native_str
 __all__ = [text_to_native_str(n) for n in __all__]
 
+@deprecate(version="3.4.3",
+           message="Pass a ``residualFn`` to :meth:`~fipy.terms.term.Term.sweep`.")
 def residual(var, matrix, RHSvector):
     r"""
     Determines the residual for the current solution matrix and variable.
@@ -34,6 +37,9 @@ def residual(var, matrix, RHSvector):
     Lx = matrix * array(var)
     return LINFnorm(Lx - RHSvector)
 
+@deprecate(version="3.4.3",
+           message="Calculate error explicitly with "
+                   ":mod:`~fipy.tools.numerix`.")
 def error(var, matrix, RHSvector, norm):
     r"""
     .. math::
@@ -60,6 +66,9 @@ def error(var, matrix, RHSvector, norm):
     denom = L1norm(var.old)
     return L1norm(var - var.old) / (denom + (denom == 0))
 
+@deprecate(version="3.4.3",
+           message="Calculate error explicitly with "
+                   ":mod:`~fipy.tools.numerix.L1norm`.")
 def L1error(var, matrix, RHSvector):
     r"""
     .. math::
@@ -81,6 +90,9 @@ def L1error(var, matrix, RHSvector):
     from fipy.tools.numerix import L1norm
     return error(var, matrix, RHSvector, L1norm)
 
+@deprecate(version="3.4.3",
+           message="Calculate error explicitly with "
+                   ":mod:`~fipy.tools.numerix.L2norm`.")
 def L2error(var, matrix, RHSvector):
     r"""
     .. math::
@@ -102,6 +114,9 @@ def L2error(var, matrix, RHSvector):
     from fipy.tools.numerix import L2norm
     return error(var, matrix, RHSvector, L2norm)
 
+@deprecate(version="3.4.3",
+           message="Calculate error explicitly with "
+                   ":mod:`~fipy.tools.numerix.LINFnorm`.")
 def LINFerror(var, matrix, RHSvector):
     r"""
     .. math::
@@ -123,6 +138,8 @@ def LINFerror(var, matrix, RHSvector):
     from fipy.tools.numerix import LINFnorm
     return error(var, matrix, RHSvector, LINFnorm)
 
+@deprecate(version="3.4.3",
+           message="Use the :term:`steppyngstounes` package instead.")
 def sweepMonotonic(fn, *args, **kwargs):
     """
     Repeatedly calls :func:`fn(*args, **kwargs)` until the residual returned by
diff --git a/fipy/steppers/pidStepper.py b/fipy/steppers/pidStepper.py
index 6de962f224..3b1ec30bb2 100644
--- a/fipy/steppers/pidStepper.py
+++ b/fipy/steppers/pidStepper.py
@@ -1,11 +1,14 @@
 from __future__ import division
 from __future__ import unicode_literals
 from fipy.steppers.stepper import Stepper
+from fipy.tools.decorators import deprecate
 
 __all__ = ["PIDStepper"]
 from future.utils import text_to_native_str
 __all__ = [text_to_native_str(n) for n in __all__]
 
+@deprecate(version="3.4.3",
+           message="Use the :term:`steppyngstounes` package instead.")
 class PIDStepper(Stepper):
     """
     Adaptive stepper using a PID controller, based on::
diff --git a/fipy/steppers/pseudoRKQSStepper.py b/fipy/steppers/pseudoRKQSStepper.py
index 20a9dd6694..8ba0d764e7 100644
--- a/fipy/steppers/pseudoRKQSStepper.py
+++ b/fipy/steppers/pseudoRKQSStepper.py
@@ -1,10 +1,13 @@
 from __future__ import unicode_literals
 from fipy.steppers.stepper import Stepper
+from fipy.tools.decorators import deprecate
 
 __all__ = ["PseudoRKQSStepper"]
 from future.utils import text_to_native_str
 __all__ = [text_to_native_str(n) for n in __all__]
 
+@deprecate(version="3.4.3",
+           message="Use the :term:`steppyngstounes` package instead.")
 class PseudoRKQSStepper(Stepper):
     """
     Adaptive stepper based on the ``rkqs`` (Runge-Kutta
diff --git a/fipy/steppers/stepper.py b/fipy/steppers/stepper.py
index dce9b3c062..4cdd2e2544 100644
--- a/fipy/steppers/stepper.py
+++ b/fipy/steppers/stepper.py
@@ -2,10 +2,14 @@
 from builtins import object
 __docformat__ = 'restructuredtext'
 
+# from fipy.tools.decorators import deprecate
+
 __all__ = ["Stepper"]
 from future.utils import text_to_native_str
 __all__ = [text_to_native_str(n) for n in __all__]
 
+# @deprecate(version="3.4.3",
+#            message="Use the :term:`steppyngstounes` package instead.")
 class Stepper(object):
     def __init__(self, vardata=()):
         self.vardata = vardata