diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py
index db906b3c13..961afefb6c 100644
--- a/openmc/deplete/integrator/__init__.py
+++ b/openmc/deplete/integrator/__init__.py
@@ -5,6 +5,7 @@ Integrator
 The integrator subcomponents.
 """
 
+from .abc import Integrator
 from .cf4 import *
 from .cecm import *
 from .celi import *
diff --git a/openmc/deplete/integrator/abc.py b/openmc/deplete/integrator/abc.py
new file mode 100644
index 0000000000..d6f318fa4c
--- /dev/null
+++ b/openmc/deplete/integrator/abc.py
@@ -0,0 +1,152 @@
+from copy import deepcopy
+from abc import ABC, abstractmethod
+from collections.abc import Iterable
+
+from openmc.deplete import Results
+
+
+class Integrator(ABC):
+    """Abstract class for solving the time-integration for depletion
+
+    """
+
+    def __init__(self, operator, timesteps, power=None, power_density=None):
+        """
+        Parameters
+        ----------
+        operator : openmc.deplete.TransportOperator
+            The operator object to simulate on.
+        timesteps : iterable of float
+            Array of timesteps in units of [s]. Note that values are not
+            cumulative.
+        power : float or iterable of float, optional
+            Power of the reactor in [W]. A single value indicates that
+            the power is constant over all timesteps. An iterable
+            indicates potentially different power levels for each timestep.
+            For a 2D problem, the power can be given in [W/cm] as long
+            as the "volume" assigned to a depletion material is actually
+            an area in [cm^2]. Either ``power`` or ``power_density`` must be
+            specified.
+        power_density : float or iterable of float, optional
+            Power density of the reactor in [W/gHM]. It is multiplied by
+            initial heavy metal inventory to get total power if ``power``
+            is not speficied.
+        """
+        self.operator = operator
+        self.chain = operator.chain
+        if not isinstance(timesteps, Iterable):
+            self.timesteps = [timesteps]
+        else:
+            self.timesteps = timesteps
+        if power is None:
+            if power_density is None:
+                raise ValueError("Either power or power density must be set")
+            if not isinstance(power_density, Iterable):
+                power = power_density * operator.heavy_metal
+            else:
+                power = [p * operator.heavy_metal for p in power_density]
+
+        if not isinstance(power, Iterable):
+            # TODO Maybe use itertools.zip_longest?
+            # Ensure that power is single value if that is the case
+            power = [power] * len(self.timesteps)
+
+        self.power = power
+
+    @abstractmethod
+    def __call__(self, conc, rates, dt, power):
+        """Perform the integration across one time step
+
+        Parameters
+        ----------
+        conc : numpy.ndarray
+            Initial concentrations for all nuclides in [atom]
+        rates : openmc.deplete.ReactionRates
+            Reaction rates from operator
+        dt : float
+            Time in [s] for the entire depletion interval
+        power : float
+            Power of the system [W]
+
+        Returns
+        -------
+        proc_time : float
+            Time spent in CRAM routines for all materials
+        conc_list : list of numpy.ndarray
+            Concentrations at each of the intermediate points with
+            the final concentration as the last element
+        op_results : list of openmc.deplete.OperatorResult
+            Eigenvalue and reaction rates from intermediate transport
+            simulations
+        """
+
+    def __iter__(self):
+        for dt, p in zip(self.timesteps, self.power):
+            yield dt, p
+
+    def __len__(self):
+        return len(self.timesteps)
+
+    def _get_bos_data(self, step_index, step_power, prev_conc):
+        if step_index > 0 or self.operator.prev_res is None:
+            x = deepcopy(prev_conc)
+            res = self.operator(x, step_power)
+        else:
+            # Get previous concentration
+            x = self.operator.prev_res[-1].data[0]
+
+            # Get reaction rates and keff
+            res = self.operator.prev_res[-1]
+            res.rates = res.rates[0]
+            res.k = res.k[0]
+
+            # Scale rates by ratio of powers
+            res.rates *= step_power / res.power[0]
+        return x, res
+
+    def _get_start_data(self):
+        if self.operator.prev_res is None:
+            return 0.0, 0
+        return self.operator.prev_res[-1].time[-1], len(self.operator.prev_res)
+
+    def integrate_all(self):
+        """Perform the entire depletion process across all steps"""
+        with self.operator as conc:
+            t, i_start = self._get_start_data()
+
+            for i, (dt, p) in enumerate(self):
+                conc, res = self._get_bos_data(i, p, conc)
+                proc_time, conc_list, res_list = self(conc, res.rates, dt, p)
+
+                # Insert BOS concentration, transport results
+                conc_list.insert(0, conc)
+                res_list.insert(0, res)
+
+                # Remove actual EOS concentration for next step
+                conc = conc_list.pop()
+
+                self._save_results(
+                    conc_list, res_list, [t, t + dt], p, i_start + i,
+                    proc_time)
+
+                t += dt
+
+            # Final simulation
+            res_list = [self.operator(conc, p)]
+            self._save_results(
+                [conc], res_list, [t, t], p, i_start + len(self))
+
+    def _save_results(self, conc_list, results_list, time_list, power,
+                      index, proc_time=None):
+        """Save the results at the end of of one step
+
+        Abstracted to support the predictor's unique save location
+        """
+        Results.save(
+            self.operator, conc_list, results_list, time_list,
+            power, index, proc_time)
+
+    @classmethod
+    def integrate(cls, operator, timesteps, power=None, power_density=None):
+        """High-level interface for depleting with this integrator"""
+        return cls(operator, timesteps, power, power_density).integrate_all()
diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py
index 758302b027..031e4ce7b1 100644
--- a/openmc/deplete/integrator/cecm.py
+++ b/openmc/deplete/integrator/cecm.py
@@ -1,13 +1,12 @@
 """The CE/CM integrator."""
 
-import copy
-from collections.abc import Iterable
+from textwrap import dedent
 
+from .abc import Integrator
 from .cram import timed_deplete
-from ..results import Results
 
 
-def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
+class CECMIntegrator(Integrator):
     r"""Deplete using the CE/CM algorithm.
 
     Implements the second order `CE/CM predictor-corrector algorithm
@@ -24,100 +23,51 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
         A_c &= A(y_m, t_n + h/2) \\
         y_{n+1} &= \text{expm}(A_c h) y_n
         \end{aligned}
-
-    Parameters
-    ----------
-    operator : openmc.deplete.TransportOperator
-        The operator object to simulate on.
-    timesteps : iterable of float
-        Array of timesteps in units of [s]. Note that values are not cumulative.
-    power : float or iterable of float, optional
-        Power of the reactor in [W]. A single value indicates that the power is
-        constant over all timesteps. An iterable indicates potentially different
-        power levels for each timestep. For a 2D problem, the power can be given
-        in [W/cm] as long as the "volume" assigned to a depletion material is
-        actually an area in [cm^2]. Either `power` or `power_density` must be
-        specified.
-    power_density : float or iterable of float, optional
-        Power density of the reactor in [W/gHM]. It is multiplied by initial
-        heavy metal inventory to get total power if `power` is not speficied.
-    print_out : bool, optional
-        Whether or not to print out time.
-
     """
-    if power is None:
-        if power_density is None:
-            raise ValueError(
-                "Neither power nor power density was specified.")
-        if not isinstance(power_density, Iterable):
-            power = power_density*operator.heavy_metal
-        else:
-            power = [i*operator.heavy_metal for i in power_density]
 
-    if not isinstance(power, Iterable):
-        power = [power]*len(timesteps)
+    def __call__(self, conc, rates, dt, power):
+        """Integrate using CE/CM
 
-    # Generate initial conditions
-    with operator as vec:
-        # Initialize time and starting index
-        if operator.prev_res is None:
-            t = 0.0
-            i_res = 0
-        else:
-            t = operator.prev_res[-1].time[-1]
-            i_res = len(operator.prev_res)
+        Parameters
+        ----------
+        conc : numpy.ndarray
+            Initial concentrations for all nuclides in [atom]
+        rates : openmc.deplete.ReactionRates
+            Reaction rates from operator
+        dt : float
+            Time in [s] for the entire depletion interval
+        power : float
+            Power of the system [W]
 
-        chain = operator.chain
+        Returns
+        -------
+        proc_time : float
+            Time spent in CRAM routines for all materials
+        conc_list : list of numpy.ndarray
+            Concentrations at each of the intermediate points with
+            the final concentration as the last element
+        op_results : list of openmc.deplete.OperatorResult
+            Eigenvalue and reaction rates from transport simulations
+        """
+        # deplete across first half of inteval
+        time0, x_middle = timed_deplete(self.chain, conc, rates, dt / 2)
+        res_middle = self.operator(x_middle, power)
 
-        for i, (dt, p) in enumerate(zip(timesteps, power)):
-            # Get beginning-of-timestep concentrations and reaction rates
-            # Avoid doing first transport run if already done in previous
-            # calculation
-            if i > 0 or operator.prev_res is None:
-                x = [copy.deepcopy(vec)]
-                op_results = [operator(x[0], p)]
+        # deplete across entire interval with BOS concentrations,
+        # MOS reaction rates
+        time1, x_end = timed_deplete(self.chain, conc, res_middle.rates, dt)
 
-            else:
-                # Get initial concentration
-                x = [operator.prev_res[-1].data[0]]
+        return time0 + time1, [x_middle, x_end], [res_middle]
 
-                # Get rates
-                op_results = [operator.prev_res[-1]]
-                op_results[0].rates = op_results[0].rates[0]
 
-                # Set first stage value of keff
-                op_results[0].k = op_results[0].k[0]
+def cecm(operator, timesteps, power=None, power_density=None, print_out=False):
+    # TODO Remove print_out since depletion timings are stored
+    return CECMIntegrator(
+        operator, timesteps, power, power_density).integrate_all()
 
-                # Scale reaction rates by ratio of powers
-                power_res = operator.prev_res[-1].power
-                ratio_power = p / power_res
-                op_results[0].rates *= ratio_power[0]
 
-            # Deplete for first half of timestep
-            proc_time, x_middle = timed_deplete(
-                chain, x[0], op_results[0].rates, dt/2, print_out)
-
-            # Get middle-of-timestep reaction rates
-            x.append(x_middle)
-            op_results.append(operator(x_middle, p))
-
-            # Deplete for full timestep using beginning-of-step materials
-            # and middle-of-timestep reaction rates
-            pt_end, x_end = timed_deplete(
-                chain, x[0], op_results[1].rates, dt, print_out)
-
-            # Create results, write to disk
-            Results.save(
-                operator, x, op_results, [t, t + dt], p, i_res + i,
-                proc_time + pt_end)
-
-            # Advance time, update vector
-            t += dt
-            vec = copy.deepcopy(x_end)
-
-        # Perform one last simulation
-        x = [copy.deepcopy(vec)]
-        op_results = [operator(x[0], power[-1])]
-
-        # Create results, write to disk
-        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
+try:
+    cecm.__doc__ = (
+        dedent(CECMIntegrator.__doc__) + dedent(Integrator.__init__.__doc__))
+except AttributeError:
+    pass