diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index d4055f0fd..532e06655 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -6,8 +6,10 @@
 
 .. module:: openmc.deplete
 
-Two functions are provided that implement different time-integration algorithms
-for depletion calculations.
+Several functions are provided that implement different time-integration
+algorithms for depletion calculations, which are described in detail in Colin
+Josey's thesis, `Development and analysis of high order neutron
+transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
 
 .. autosummary::
    :toctree: generated
@@ -16,6 +18,12 @@ for depletion calculations.
 
    integrator.predictor
    integrator.cecm
+   integrator.celi
+   integrator.leqi
+   integrator.cf4
+   integrator.epc_rk4
+   integrator.si_celi
+   integrator.si_leqi
 
 Each of these functions expects a "transport operator" to be passed. An operator
 specific to OpenMC is available using the following class:
diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py
index cf8caffdf..db906b3c1 100644
--- a/openmc/deplete/integrator/__init__.py
+++ b/openmc/deplete/integrator/__init__.py
@@ -5,6 +5,12 @@ Integrator
 The integrator subcomponents.
 """
 
+from .cf4 import *
 from .cecm import *
+from .celi import *
 from .cram import *
+from .epc_rk4 import *
+from .leqi import *
 from .predictor import *
+from .si_celi import *
+from .si_leqi import *
diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py
index ac519cf87..42024863a 100644
--- a/openmc/deplete/integrator/cecm.py
+++ b/openmc/deplete/integrator/cecm.py
@@ -11,8 +11,10 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
     r"""Deplete using the CE/CM algorithm.
 
     Implements the second order `CE/CM predictor-corrector algorithm
-    <https://doi.org/10.13182/NSE14-92>`_.  This algorithm is mathematically
-    defined as:
+    <https://doi.org/10.13182/NSE14-92>`_.
+
+    "CE/CM" stands for constant extrapolation on predictor and constant
+    midpoint on corrector. This algorithm is mathematically defined as:
 
     .. math::
         y' &= A(y, t) y(t)
@@ -59,20 +61,16 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
 
     # Generate initial conditions
     with operator as vec:
-        chain = operator.chain
-
-        # Initialize time
+        # Initialize time and starting index
         if operator.prev_res is None:
             t = 0.0
-        else:
-            t = operator.prev_res[-1].time[-1]
-
-        # Initialize starting index for saving results
-        if operator.prev_res is None:
             i_res = 0
         else:
+            t = operator.prev_res[-1].time[-1]
             i_res = len(operator.prev_res)
 
+        chain = operator.chain
+
         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
@@ -95,10 +93,10 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
                 # Scale reaction rates by ratio of powers
                 power_res = operator.prev_res[-1].power
                 ratio_power = p / power_res
-                op_results[0].rates[0] *= ratio_power[0]
+                op_results[0].rates *= ratio_power[0]
 
             # Deplete for first half of timestep
-            x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)
+            x_middle = deplete(chain, x[0], op_results[0].rates, dt/2, print_out)
 
             # Get middle-of-timestep reaction rates
             x.append(x_middle)
@@ -106,7 +104,7 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
 
             # Deplete for full timestep using beginning-of-step materials
             # and middle-of-timestep reaction rates
-            x_end = deplete(chain, x[0], op_results[1], dt, print_out)
+            x_end = 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)
diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py
new file mode 100644
index 000000000..0ba86d585
--- /dev/null
+++ b/openmc/deplete/integrator/celi.py
@@ -0,0 +1,166 @@
+"""The CE/LI CFQ4 integrator."""
+
+import copy
+from collections.abc import Iterable
+
+from .cram import deplete
+from ..results import Results
+
+
+# Functions to form the special matrix for depletion
+def _celi_f1(chain, rates):
+    return 5/12 * chain.form_matrix(rates[0]) + \
+           1/12 * chain.form_matrix(rates[1])
+
+def _celi_f2(chain, rates):
+    return 1/12 * chain.form_matrix(rates[0]) + \
+           5/12 * chain.form_matrix(rates[1])
+
+def celi(operator, timesteps, power=None, power_density=None,
+         print_out=True):
+    r"""Deplete using the CE/LI CFQ4 algorithm.
+
+    Implements the CE/LI Predictor-Corrector algorithm using the `fourth order
+    commutator-free integrator <https://doi.org/10.1137/05063042>`_.
+
+    "CE/LI" stands for constant extrapolation on predictor and linear
+    interpolation on corrector. This algorithm is mathematically defined as:
+
+    .. math::
+        y' &= A(y, t) y(t)
+
+        A_0 &= A(y_n, t_n)
+
+        y_p &= \text{expm}(h A_0) y_n
+
+        A_1 &= A(y_p, t_n + h)
+
+        y_{n+1} &= \text{expm}(\frac{h}{12} A_0 + \frac{5h}{12} A1)
+                   \text{expm}(\frac{5h}{12} A_0 + \frac{h}{12} A1) y_n
+
+    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)
+
+    # 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)
+
+        for i, (dt, p) in enumerate(zip(timesteps, power)):
+            vec, t, _ = celi_inner(operator, vec, p, i, i_res, t, dt,
+                                   print_out)
+
+        # 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))
+
+def celi_inner(operator, vec, p, i, i_res, t, dt, print_out):
+    """ The inner loop of CE/LI CFQ4.
+
+    Parameters
+    ----------
+    operator : Operator
+        The operator object to simulate on.
+    x : list of nuclide vector
+        Nuclide vector, beginning of time.
+    p : float
+        Power of the reactor in [W]
+    i : int
+        Current iteration number.
+    i_res : int
+        Starting index, for restart calculation.
+    t : float
+        Time at start of step.
+    dt : float
+        Time step.
+    print_out : bool
+        Whether or not to print out time.
+
+    Returns
+    -------
+    list of numpy.array
+        Nuclide vector, end of time.
+    float
+        Next time
+    OperatorResult
+        Operator result from beginning of step.
+    """
+
+    chain = operator.chain
+
+    # 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)]
+
+    else:
+        # Get initial concentration
+        x = [operator.prev_res[-1].data[0]]
+
+        # 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]
+
+        # 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 to end
+    x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out)
+    x.append(x_new)
+    op_results.append(operator(x[1], p))
+
+    # Deplete with two matrix exponentials
+    rates = list(zip(op_results[0].rates, op_results[1].rates))
+    x_end = deplete(chain, x[0], rates, dt, print_out,
+                    matrix_func=_celi_f1)
+    x_end = deplete(chain, x_end, rates, dt, print_out,
+                    matrix_func=_celi_f2)
+
+    # Create results, write to disk
+    Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+
+    # return updated time and vectors
+    return x_end, t + dt, op_results[0]
diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py
new file mode 100644
index 000000000..e93b22d5c
--- /dev/null
+++ b/openmc/deplete/integrator/cf4.py
@@ -0,0 +1,161 @@
+"""The CF4 integrator."""
+
+import copy
+from collections.abc import Iterable
+
+from .cram import deplete
+from ..results import Results
+
+
+# Functions to form the special matrix for depletion
+def _cf4_f1(chain, rates):
+    return 1/2 * chain.form_matrix(rates)
+
+def _cf4_f2(chain, rates):
+    return -1/2 * chain.form_matrix(rates[0]) + \
+                  chain.form_matrix(rates[1])
+
+def _cf4_f3(chain, rates):
+    return  1/4  * chain.form_matrix(rates[0]) + \
+            1/6  * chain.form_matrix(rates[1]) + \
+            1/6  * chain.form_matrix(rates[2]) + \
+           -1/12 * chain.form_matrix(rates[3])
+
+def _cf4_f4(chain, rates):
+    return -1/12 * chain.form_matrix(rates[0]) + \
+            1/6  * chain.form_matrix(rates[1]) + \
+            1/6  * chain.form_matrix(rates[2]) + \
+            1/4  * chain.form_matrix(rates[3])
+
+def cf4(operator, timesteps, power=None, power_density=None, print_out=True):
+    r"""Deplete using the CF4 algorithm.
+
+    Implements the fourth order `commutator-free Lie algorithm
+    <https://doi.org/10.1016/S0167-739X(02)00161-9>`_.
+    This algorithm is mathematically defined as:
+
+    .. math::
+        F_1 &= h A(y_0)
+
+        y_1 &= \text{expm}(1/2 F_1) y_0
+
+        F_2 &= h A(y_1)
+
+        y_2 &= \text{expm}(1/2 F_2) y_0
+
+        F_3 &= h A(y_2)
+
+        y_3 &= \text{expm}(-1/2 F_1 + F_3) y_1
+
+        F_4 &= h A(y_3)
+
+        y_4 &= \text{expm}( 1/4  F_1 + 1/6 F_2 + 1/6 F_3 - 1/12 F_4)
+               \text{expm}(-1/12 F_1 + 1/6 F_2 + 1/6 F_3 + 1/4  F_4) y_0
+
+    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)
+
+    # 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)
+
+        chain = operator.chain
+
+        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)]
+
+            else:
+                # Get initial concentration
+                x = [operator.prev_res[-1].data[0]]
+
+                # 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]
+
+                # 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]
+
+            # Step 1: deplete with matrix 1/2*A(y0)
+            x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out,
+                            matrix_func=_cf4_f1)
+            x.append(x_new)
+            op_results.append(operator(x_new, p))
+
+            # Step 2: deplete with matrix 1/2*A(y1)
+            x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out,
+                            matrix_func=_cf4_f1)
+            x.append(x_new)
+            op_results.append(operator(x_new, p))
+
+            # Step 3: deplete with matrix -1/2*A(y0)+A(y2)
+            rates = list(zip(op_results[0].rates, op_results[2].rates))
+            x_new = deplete(chain, x[1], rates, dt, print_out,
+                            matrix_func=_cf4_f2)
+            x.append(x_new)
+            op_results.append(operator(x_new, p))
+
+            # Step 4: deplete with two matrix exponentials
+            rates = list(zip(op_results[0].rates, op_results[1].rates,
+                             op_results[2].rates, op_results[3].rates))
+            x_end = deplete(chain, x[0], rates, dt, print_out,
+                            matrix_func=_cf4_f3)
+            x_end = deplete(chain, x_end, rates, dt, print_out,
+                            matrix_func=_cf4_f4)
+
+            # Create results, write to disk
+            Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+
+            # 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))
diff --git a/openmc/deplete/integrator/cram.py b/openmc/deplete/integrator/cram.py
index 85954603a..8a8940df8 100644
--- a/openmc/deplete/integrator/cram.py
+++ b/openmc/deplete/integrator/cram.py
@@ -14,7 +14,7 @@ import scipy.sparse.linalg as sla
 from .. import comm
 
 
-def deplete(chain, x, op_result, dt, print_out):
+def deplete(chain, x, rates, dt, print_out=True, matrix_func=None):
     """Deplete materials using given reaction rates for a specified time
 
     Parameters
@@ -23,12 +23,14 @@ def deplete(chain, x, op_result, dt, print_out):
         Depletion chain
     x : list of numpy.ndarray
         Atom number vectors for each material
-    op_result : openmc.deplete.OperatorResult
-        Result of applying transport operator (contains reaction rates)
+    rates : openmc.deplete.ReactionRates
+        Reaction rates (from transport operator)
     dt : float
         Time in [s] to deplete for
-    print_out : bool
+    print_out : bool, optional
         Whether to show elapsed time
+    maxtrix_func : function, optional
+        Function to form the depletion matrix
 
     Returns
     -------
@@ -38,16 +40,9 @@ def deplete(chain, x, op_result, dt, print_out):
     """
     t_start = time.time()
 
-    # Set up iterators
-    n_mats = len(x)
-    chains = repeat(chain, n_mats)
-    vecs = (x[i] for i in range(n_mats))
-    rates = (op_result.rates[i, :, :] for i in range(n_mats))
-    dts = repeat(dt, n_mats)
-
     # Use multiprocessing pool to distribute work
     with Pool() as pool:
-        iters = zip(chains, vecs, rates, dts)
+        iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func))
         x_result = list(pool.starmap(_cram_wrapper, iters))
 
     t_end = time.time()
@@ -58,12 +53,12 @@ def deplete(chain, x, op_result, dt, print_out):
     return x_result
 
 
-def _cram_wrapper(chain, n0, rates, dt):
+def _cram_wrapper(chain, n0, rates, dt, matrix_func=None):
     """Wraps depletion matrix creation / CRAM solve for multiprocess execution
 
     Parameters
     ----------
-    chain : DepletionChain
+    chain : openmc.deplete.Chain
         Depletion chain used to construct the burnup matrix
     n0 : numpy.array
         Vector to operate a matrix exponent on.
@@ -71,13 +66,19 @@ def _cram_wrapper(chain, n0, rates, dt):
         2D array indexed by nuclide then by cell.
     dt : float
         Time to integrate to.
+    maxtrix_func : function, optional
+        Function to form the depletion matrix
 
     Returns
     -------
     numpy.array
         Results of the matrix exponent.
     """
-    A = chain.form_matrix(rates)
+
+    if matrix_func is None:
+        A = chain.form_matrix(rates)
+    else:
+        A = matrix_func(chain, rates)
     return CRAM48(A, n0, dt)
 
 
diff --git a/openmc/deplete/integrator/epc_rk4.py b/openmc/deplete/integrator/epc_rk4.py
new file mode 100644
index 000000000..3789f3698
--- /dev/null
+++ b/openmc/deplete/integrator/epc_rk4.py
@@ -0,0 +1,147 @@
+"""The EPC-RK4 integrator."""
+
+import copy
+from collections.abc import Iterable
+
+from .cram import deplete
+from ..results import Results
+
+
+# Functions to form the special matrix for depletion
+def _rk4_f1(chain, rates):
+    return 1/2 * chain.form_matrix(rates)
+
+def _rk4_f4(chain, rates):
+    return 1/6 * chain.form_matrix(rates[0]) + \
+           1/3 * chain.form_matrix(rates[1]) + \
+           1/3 * chain.form_matrix(rates[2]) + \
+           1/6 * chain.form_matrix(rates[3])
+
+def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True):
+    r"""Deplete using the EPC-RK4 algorithm.
+
+    Implements an extended predictor-corrector algorithm with traditional
+    Runge-Kutta 4 method.
+    This algorithm is mathematically defined as:
+
+    .. math::
+        F_1 &= h A(y_0)
+
+        y_1 &= \text{expm}(1/2 F_1) y_0
+
+        F_2 &= h A(y_1)
+
+        y_2 &= \text{expm}(1/2 F_2) y_0
+
+        F_3 &= h A(y_2)
+
+        y_3 &= \text{expm}(F_3) y_0
+
+        F_4 &= h A(y_3)
+
+        y_4 &= \text{expm}(1/6 F_1 + 1/3 F_2 + 1/3 F_3 + 1/6 F_4) y_0
+
+    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)
+
+    # 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)
+
+        chain = operator.chain
+
+        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)]
+
+            else:
+                # Get initial concentration
+                x = [operator.prev_res[-1].data[0]]
+
+                # 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]
+
+                # 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]
+
+            # Step 1: deplete with matrix 1/2*A(y0)
+            x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out,
+                            matrix_func=_rk4_f1)
+            x.append(x_new)
+            op_results.append(operator(x[1], p))
+
+            # Step 2: deplete with matrix 1/2*A(y1)
+            x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out,
+                            matrix_func=_rk4_f1)
+            x.append(x_new)
+            op_results.append(operator(x[2], p))
+
+            # Step 3: deplete with matrix A(y2)
+            x_new = deplete(chain, x[0], op_results[2].rates, dt, print_out)
+            x.append(x_new)
+            op_results.append(operator(x[3], p))
+
+            # Step 4: deplete with matrix 1/6*A(y0)+1/3*A(y1)+1/3*A(y2)+1/6*A(y3)
+            rates = list(zip(op_results[0].rates, op_results[1].rates,
+                             op_results[2].rates, op_results[3].rates))
+            x_end = deplete(chain, x[0], rates, dt, print_out,
+                            matrix_func=_rk4_f4)
+
+            # Create results, write to disk
+            Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+
+            # 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))
diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py
new file mode 100644
index 000000000..9d221e877
--- /dev/null
+++ b/openmc/deplete/integrator/leqi.py
@@ -0,0 +1,170 @@
+"""The LE/QI CFQ4 integrator."""
+
+import copy
+from collections.abc import Iterable
+from itertools import repeat
+
+from .celi import celi_inner
+from .cram import deplete
+from ..results import Results
+
+
+# Functions to form the special matrix for depletion
+def _leqi_f1(chain, inputs):
+    f1 = chain.form_matrix(inputs[0])
+    f2 = chain.form_matrix(inputs[1])
+    dt_l, dt = inputs[2], inputs[3]
+    return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2
+
+def _leqi_f2(chain, inputs):
+    f1 = chain.form_matrix(inputs[0])
+    f2 = chain.form_matrix(inputs[1])
+    dt_l, dt = inputs[2], inputs[3]
+    return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2
+
+def _leqi_f3(chain, inputs):
+    f1 = chain.form_matrix(inputs[0])
+    f2 = chain.form_matrix(inputs[1])
+    f3 = chain.form_matrix(inputs[2])
+    dt_l, dt = inputs[3], inputs[4]
+    return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \
+           (dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \
+           dt_l / (12 * (dt + dt_l)) * f3
+
+def _leqi_f4(chain, inputs):
+    f1 = chain.form_matrix(inputs[0])
+    f2 = chain.form_matrix(inputs[1])
+    f3 = chain.form_matrix(inputs[2])
+    dt_l, dt = inputs[3], inputs[4]
+    return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \
+           (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \
+           (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3
+
+def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
+    r"""Deplete using the LE/QI CFQ4 algorithm.
+
+    Implements the LE/QI Predictor-Corrector algorithm using the `fourth order
+    commutator-free integrator <https://doi.org/10.1137/05063042>`_.
+
+    "LE/QI" stands for linear extrapolation on predictor and quadratic
+    interpolation on corrector. This algorithm is mathematically defined as:
+
+    .. math::
+        y' &= A(y, t) y(t)
+
+        A_{last} &= A(y_{n-1}, t_n - h_1)
+
+        A_0 &= A(y_n, t_n)
+
+        F_1 &= \frac{-h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+h_2)}{12h_1} A_0
+
+        F_2 &= \frac{-5h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+5h_2)}{12h_1} A_0
+
+        y_p &= \text{expm}(F_2) \text{expm}(F_1) y_n
+
+        A_1 &= A(y_p, t_n + h_2)
+
+        F_3 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
+              \frac{h_2 (5 h_1^2 + 6 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
+              \frac{h_2 h_1)}{12 (h_1 + h_2)} A_1
+
+        F_4 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} +
+              \frac{h_2 (h_1^2 + 2 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 +
+              \frac{h_2 (5 h_1^2 + 4 h_2 h_1)}{12 h_1 (h_1 + h_2)} A_1
+
+        y_{n+1} &= \text{expm}(F_4) \text{expm}(F_3) y_n
+
+    It is initialized using the CE/LI algorithm.
+
+    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)
+
+    # 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)
+
+        chain = operator.chain
+
+        for i, (dt, p) in enumerate(zip(timesteps, power)):
+            # LE/QI needs the last step results to start
+            # Perform CE/LI CFQ4 or restore results for the first step
+            if i == 0:
+                if i_res <= 1:
+                    dt_l = dt
+                    x_new, t, op_res_last = celi_inner(operator, vec, p, i,
+                                                       i_res, t, dt, print_out)
+                    continue
+                else:
+                    dt_l = t - operator.prev_res[-2].time[0]
+                    op_res_last = operator.prev_res[-2]
+                    op_res_last.rates = op_res_last.rates[0]
+                    x_new = operator.prev_res[-1].data[0]
+
+            # Perform remaining LE/QI
+            x = [copy.deepcopy(x_new)]
+            op_results = [operator(x[0], p)]
+
+            inputs = list(zip(op_res_last.rates, op_results[0].rates,
+                              repeat(dt_l), repeat(dt)))
+            x_new = deplete(chain, x[0], inputs, dt, print_out,
+                            matrix_func=_leqi_f1)
+            x_new = deplete(chain, x_new, inputs, dt, print_out,
+                            matrix_func=_leqi_f2)
+            x.append(x_new)
+            op_results.append(operator(x[1], p))
+
+            inputs = list(zip(op_res_last.rates, op_results[0].rates,
+                              op_results[1].rates, repeat(dt_l), repeat(dt)))
+            x_new = deplete(chain, x[0], inputs, dt, print_out,
+                            matrix_func=_leqi_f3)
+            x_new = deplete(chain, x_new, inputs, dt, print_out,
+                            matrix_func=_leqi_f4)
+
+            # Create results, write to disk
+            Results.save(operator, x, op_results, [t, t+dt], p, i_res+i)
+
+            # update results
+            op_res_last = copy.deepcopy(op_results[0])
+            t += dt
+            dt_l = dt
+
+        # Perform one last simulation
+        x = [copy.deepcopy(x_new)]
+        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))
diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py
index 969144f68..f670a125e 100644
--- a/openmc/deplete/integrator/predictor.py
+++ b/openmc/deplete/integrator/predictor.py
@@ -55,20 +55,16 @@ def predictor(operator, timesteps, power=None, power_density=None,
 
     # Generate initial conditions
     with operator as vec:
-        chain = operator.chain
-
-        # Initialize time
+        # Initialize time and starting index
         if operator.prev_res is None:
             t = 0.0
-        else:
-            t = operator.prev_res[-1].time[-1]
-
-        # Initialize starting index for saving results
-        if operator.prev_res is None:
             i_res = 0
         else:
+            t = operator.prev_res[-1].time[-1]
             i_res = len(operator.prev_res) - 1
 
+        chain = operator.chain
+
         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
@@ -90,10 +86,10 @@ def predictor(operator, timesteps, power=None, power_density=None,
                 # Scale reaction rates by ratio of powers
                 power_res = operator.prev_res[-1].power
                 ratio_power = p / power_res
-                op_results[0].rates[0] *= ratio_power[0]
+                op_results[0].rates *= ratio_power[0]
 
             # Deplete for full timestep
-            x_end = deplete(chain, x[0], op_results[0], dt, print_out)
+            x_end = deplete(chain, x[0], op_results[0].rates, dt, print_out)
 
             # Advance time, update vector
             t += dt
diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py
new file mode 100644
index 000000000..4f34bace6
--- /dev/null
+++ b/openmc/deplete/integrator/si_celi.py
@@ -0,0 +1,163 @@
+"""The SI-CE/LI CFQ4 integrator."""
+
+import copy
+from collections.abc import Iterable
+
+from .cram import deplete
+from ..results import Results
+from ..abc import OperatorResult
+from .celi import _celi_f1, _celi_f2
+
+
+def si_celi(operator, timesteps, power=None, power_density=None,
+            print_out=True, m=10):
+    r"""Deplete using the SI-CE/LI CFQ4 algorithm.
+
+    Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using
+    the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
+
+    Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
+    <http://hdl.handle.net/1721.1/113721>`_.
+
+    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.
+    m : int, optional
+        Number of stages.
+    """
+    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)
+
+    # 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)
+
+        # Get the concentrations and reaction rates for the first
+        # beginning-of-timestep (BOS). Compute with m (stage number) times as
+        # many neutrons as later simulations for statistics reasons if no
+        # previous calculation results present
+        if operator.prev_res is None:
+            x = [copy.deepcopy(vec)]
+            if hasattr(operator, "settings"):
+                operator.settings.particles *= m
+            op_results = [operator(x[0], power[0])]
+            if hasattr(operator, "settings"):
+                operator.settings.particles //= m
+        else:
+            # Get initial concentration
+            x = [operator.prev_res[-1].data[0]]
+
+            # 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]
+
+            # Scale reaction rates by ratio of powers
+            power_res = operator.prev_res[-1].power
+            ratio_power = power[0] / power_res
+            op_results[0].rates *= ratio_power[0]
+
+        for i, (dt, p) in enumerate(zip(timesteps, power)):
+            x, t, op_results = si_celi_inner(operator, x, op_results, p,
+                                             i, i_res, t, dt, print_out, m)
+
+        # Create results for last point, write to disk
+        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
+
+
+def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):
+    """ The inner loop of SI-CE/LI CFQ4.
+
+    Parameters
+    ----------
+    operator : Operator
+        The operator object to simulate on.
+    x : list of nuclide vector
+        Nuclide vector, beginning of time.
+    op_results : list of OperatorResult
+        Operator result at BOS.
+    p : float
+        Power of the reactor in [W]
+    i : int
+        Current iteration number.
+    i_res : int
+        Starting index, for restart calculation.
+    t : float
+        Time at start of step.
+    dt : float
+        Time step.
+    print_out : bool
+        Whether or not to print out time.
+    m : int, optional
+        Number of stages.
+
+    Returns
+    -------
+    list of nuclide vector (numpy.array)
+        Nuclide vector, end of time.
+    float
+        Next time
+    list of OperatorResult
+        Operator result at end of time.
+    """
+
+    chain = operator.chain
+
+    # Deplete to end
+    x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out)
+    x.append(x_new)
+
+    for j in range(m + 1):
+        op_res = operator(x_new, p)
+
+        if j <= 1:
+            op_res_bar = copy.deepcopy(op_res)
+        else:
+            rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates
+            k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k
+            op_res_bar = OperatorResult(k, rates)
+
+        rates = list(zip(op_results[0].rates, op_res_bar.rates))
+        x_new = deplete(chain, x[0], rates, dt, print_out,
+                        matrix_func=_celi_f1)
+        x_new = deplete(chain, x_new, rates, dt, print_out,
+                        matrix_func=_celi_f2)
+
+    # Create results, write to disk
+    op_results.append(op_res_bar)
+    Results.save(operator, x, op_results, [t, t+dt], p, i_res+i)
+
+    # return updated time and vectors
+    return [x_new], t + dt, [op_res_bar]
diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py
new file mode 100644
index 000000000..05002880a
--- /dev/null
+++ b/openmc/deplete/integrator/si_leqi.py
@@ -0,0 +1,151 @@
+"""The SI-LE/QI CFQ4 integrator."""
+
+import copy
+from collections.abc import Iterable
+from itertools import repeat
+
+from .si_celi import si_celi_inner
+from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4
+from .cram import deplete
+from ..results import Results
+from ..abc import OperatorResult
+
+
+def si_leqi(operator, timesteps, power=None, power_density=None,
+            print_out=True, m=10):
+    r"""Deplete using the SI-LE/QI CFQ4 algorithm.
+
+    Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using
+    the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
+
+    Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
+    <http://hdl.handle.net/1721.1/113721>`_.
+
+    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.
+    m : int, optional
+        Number of stages.
+    """
+    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)
+
+    # 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)
+
+        # Get the concentrations and reaction rates for the first
+        # beginning-of-timestep (BOS). Compute with m (stage number) times as
+        # many neutrons as later simulations for statistics reasons if no
+        # previous calculation results present
+        if operator.prev_res is None:
+            x = [copy.deepcopy(vec)]
+            if hasattr(operator, "settings"):
+                operator.settings.particles *= m
+            op_results = [operator(x[0], power[0])]
+            if hasattr(operator, "settings"):
+                operator.settings.particles //= m
+        else:
+            # Get initial concentration
+            x = [operator.prev_res[-1].data[0]]
+
+            # 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]
+
+            # Scale reaction rates by ratio of powers
+            power_res = operator.prev_res[-1].power
+            ratio_power = power[0] / power_res
+            op_results[0].rates *= ratio_power[0]
+
+        chain = operator.chain
+
+        for i, (dt, p) in enumerate(zip(timesteps, power)):
+            # LE/QI needs the last step results to start
+            # Perform SI-CE/LI CFQ4 or restore results for the first step
+            if i == 0:
+                dt_l = dt
+                if i_res <= 1:
+                    op_res_last = copy.deepcopy(op_results[0])
+                    x, t, op_results = si_celi_inner(operator, x, op_results, p,
+                                                     i, i_res, t, dt, print_out)
+                    continue
+                else:
+                    dt_l = t - operator.prev_res[-2].time[0]
+                    op_res_last = operator.prev_res[-2]
+                    op_res_last.rates = op_res_last.rates[0]
+                    x = [operator.prev_res[-1].data[0]]
+
+            # Perform remaining LE/QI
+            inputs = list(zip(op_res_last.rates, op_results[0].rates,
+                              repeat(dt_l), repeat(dt)))
+            x_new = deplete(chain, x[0], inputs, dt, print_out,
+                            matrix_func=_leqi_f1)
+            x_new = deplete(chain, x_new, inputs, dt, print_out,
+                            matrix_func=_leqi_f2)
+            x.append(x_new)
+
+            # Loop on inner
+            for j in range(m + 1):
+                op_res = operator(x_new, p)
+
+                if j <= 1:
+                    op_res_bar = copy.deepcopy(op_res)
+                else:
+                    rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates
+                    k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k
+                    op_res_bar = OperatorResult(k, rates)
+
+                inputs = list(zip(op_res_last.rates, op_results[0].rates,
+                                  op_res_bar.rates, repeat(dt_l), repeat(dt)))
+                x_new = deplete(chain, x[0], inputs, dt, print_out,
+                                matrix_func=_leqi_f3)
+                x_new = deplete(chain, x_new, inputs, dt, print_out,
+                                matrix_func=_leqi_f4)
+
+            # Create results, write to disk
+            op_results.append(op_res_bar)
+            Results.save(operator, x, op_results, [t, t+dt], p, i_res+i)
+
+            # update results
+            x = [x_new]
+            op_res_last = copy.deepcopy(op_results[0])
+            op_results = [op_res_bar]
+            t += dt
+            dt_l = dt
+
+        # Create results for last point, write to disk
+        Results.save(operator, x, op_results, [t, t], p, i_res+len(timesteps))
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 754e22f7d..66cac392d 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -125,6 +125,10 @@ class Operator(TransportOperator):
         else:
             self.prev_res = None
 
+        # Differentiate burnable materials with multiple instances
+        if self.diff_burnable_mats:
+            self._differentiate_burnable_mats()
+
         # Clear out OpenMC, create task lists, distribute
         openmc.reset_auto_ids()
         self.burnable_mats, volume, nuclides = self._get_burnable_mats()
@@ -227,10 +231,6 @@ class Operator(TransportOperator):
 
         """
 
-        if self.diff_burnable_mats:
-            # Automatically distribute burnable materials
-            self._differentiate_burnable_mats()
-
         burnable_mats = set()
         model_nuclides = set()
         volume = OrderedDict()
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 28d19713c..b6a89792d 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -122,7 +122,7 @@ class Model(object):
             for plot in plots:
                 self._plots.append(plot)
 
-    def deplete(self, timesteps, power, chain_file=None, method='cecm',
+    def deplete(self, timesteps, chain_file=None, method='cecm',
                 **kwargs):
         """Deplete model using specified timesteps/power
 
@@ -131,17 +131,11 @@ class Model(object):
         timesteps : iterable of float
             Array of timesteps in units of [s]. Note that values are not
             cumulative.
-        power : float or iterable of float
-            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].
         chain_file : str, optional
             Path to the depletion chain XML file.  Defaults to the
             :envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
-        method : {'cecm', 'predictor'}
-             Integration method used for depletion
+        method : str
+             Integration method used for depletion (e.g., 'cecm', 'predictor')
         **kwargs
             Keyword arguments passed to integration function (e.g.,
             :func:`openmc.deplete.integrator.cecm`)
@@ -156,12 +150,9 @@ class Model(object):
         op = dep.Operator(self.geometry, self.settings, chain_file)
 
         # Perform depletion
-        if method == 'predictor':
-            dep.integrator.predictor(op, timesteps, power, **kwargs)
-        elif method == 'cecm':
-            dep.integrator.cecm(op, timesteps, power, **kwargs)
-        else:
-            check_value('method', method, ('cecm', 'predictor'))
+        check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4',
+                                       'si_celi', 'si_leqi', 'celi', 'leqi'))
+        getattr(dep.integrator, method)(op, timesteps, **kwargs)
 
     def export_to_xml(self):
         """Export model to XML files."""
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index db7343172..b64a1681a 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -905,6 +905,7 @@ extern "C" int openmc_load_nuclide(const char* name)
       // Get filename for library containing nuclide
       int idx = it->second;
       std::string& filename = data::libraries[idx].path_;
+      write_message("Reading " + std::string{name} + " from " + filename, 6);
 
       // Open file and make sure version is sufficient
       hid_t file_id = file_open(filename, 'r');
@@ -924,6 +925,9 @@ extern "C" int openmc_load_nuclide(const char* name)
 
       // Initialize nuclide grid
       data::nuclides.back()->init_grid();
+
+      // Read multipole file into the appropriate entry on the nuclides array
+      if (settings::temperature_multipole) read_multipole_data(i_nuclide);
     } else {
       set_errmsg("Nuclide '" + std::string{name} + "' is not present in library.");
       return OPENMC_E_DATA;
diff --git a/tests/unit_tests/test_deplete_celi.py b/tests/unit_tests/test_deplete_celi.py
new file mode 100644
index 000000000..0b23b1a43
--- /dev/null
+++ b/tests/unit_tests/test_deplete_celi.py
@@ -0,0 +1,37 @@
+"""Regression tests for openmc.deplete.integrator.celi algorithm.
+
+These tests integrate a simple test problem described in dummy_geometry.py.
+"""
+
+from pytest import approx
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+def test_celi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using celi"""
+
+    op = dummy_operator.DummyOperator()
+    op.output_dir = "test_integrator_regression"
+
+    # Perform simulation using the celi algorithm
+    dt = [0.75, 0.75]
+    power = 1.0
+    openmc.deplete.celi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [1.82078767, 0.97122898]
+    s2 = [2.68441779, 0.05125966]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[2] == approx(s2[0])
+    assert y2[2] == approx(s2[1])
diff --git a/tests/unit_tests/test_deplete_cf4.py b/tests/unit_tests/test_deplete_cf4.py
new file mode 100644
index 000000000..784d228be
--- /dev/null
+++ b/tests/unit_tests/test_deplete_cf4.py
@@ -0,0 +1,37 @@
+"""Regression tests for openmc.deplete.integrator.cf4 algorithm.
+
+These tests integrate a simple test problem described in dummy_geometry.py.
+"""
+
+from pytest import approx
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+def test_cf4(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using CF4"""
+
+    op = dummy_operator.DummyOperator()
+    op.output_dir = "test_integrator_regression"
+
+    # Perform simulation using the cf4 algorithm
+    dt = [0.75, 0.75]
+    power = 1.0
+    openmc.deplete.cf4(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.06101629, 1.37783588]
+    s2 = [2.57241318, 2.63731630]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[2] == approx(s2[0])
+    assert y2[2] == approx(s2[1])
diff --git a/tests/unit_tests/test_deplete_epc_rk4.py b/tests/unit_tests/test_deplete_epc_rk4.py
new file mode 100644
index 000000000..dd15bf8e1
--- /dev/null
+++ b/tests/unit_tests/test_deplete_epc_rk4.py
@@ -0,0 +1,37 @@
+"""Regression tests for openmc.deplete.integrator.epc_rk4 algorithm.
+
+These tests integrate a simple test problem described in dummy_geometry.py.
+"""
+
+from pytest import approx
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+def test_epc_rk4(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using epc_rk4"""
+
+    op = dummy_operator.DummyOperator()
+    op.output_dir = "test_integrator_regression"
+
+    # Perform simulation using the epc_rk4 algorithm
+    dt = [0.75, 0.75]
+    power = 1.0
+    openmc.deplete.epc_rk4(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.01978516, 1.42038037]
+    s2 = [2.05246421, 3.06177191]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[2] == approx(s2[0])
+    assert y2[2] == approx(s2[1])
diff --git a/tests/unit_tests/test_deplete_leqi.py b/tests/unit_tests/test_deplete_leqi.py
new file mode 100644
index 000000000..a1a32e092
--- /dev/null
+++ b/tests/unit_tests/test_deplete_leqi.py
@@ -0,0 +1,37 @@
+"""Regression tests for openmc.deplete.integrator.leqi algorithm.
+
+These tests integrate a simple test problem described in dummy_geometry.py.
+"""
+
+from pytest import approx
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+def test_leqi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using leqi"""
+
+    op = dummy_operator.DummyOperator()
+    op.output_dir = "test_integrator_regression"
+
+    # Perform simulation using the leqi algorithm
+    dt = [0.75, 0.75]
+    power = 1.0
+    openmc.deplete.leqi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [1.82078767, 0.97122898]
+    s2 = [2.74526197, 0.23339915]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[2] == approx(s2[0])
+    assert y2[2] == approx(s2[1])
diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py
index 745f2cce5..ec62064a4 100644
--- a/tests/unit_tests/test_deplete_restart.py
+++ b/tests/unit_tests/test_deplete_restart.py
@@ -166,3 +166,236 @@ def test_restart_cecm_predictor(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+
+def test_restart_cf4(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using CF4."""
+
+    op = dummy_operator.DummyOperator()
+    output_dir = "test_restart_cf4"
+    op.output_dir = output_dir
+
+    # Perform simulation
+    dt = [0.75]
+    power = 1.0
+    openmc.deplete.cf4(op, dt, power, print_out=False)
+
+    # Load the files
+    prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    # Re-create depletion operator and load previous results
+    op = dummy_operator.DummyOperator(prev_res)
+    op.output_dir = output_dir
+
+    # Perform restarts simulation
+    openmc.deplete.cf4(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.06101629, 1.37783588]
+    s2 = [2.57241318, 2.63731630]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[3] == approx(s2[0])
+    assert y2[3] == approx(s2[1])
+
+
+def test_restart_epc_rk4(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using EPC-RK4."""
+
+    op = dummy_operator.DummyOperator()
+    output_dir = "test_restart_epc_rk4"
+    op.output_dir = output_dir
+
+    # Perform simulation
+    dt = [0.75]
+    power = 1.0
+    openmc.deplete.epc_rk4(op, dt, power, print_out=False)
+
+    # Load the files
+    prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    # Re-create depletion operator and load previous results
+    op = dummy_operator.DummyOperator(prev_res)
+    op.output_dir = output_dir
+
+    # Perform restarts simulation
+    openmc.deplete.epc_rk4(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.01978516, 1.42038037]
+    s2 = [2.05246421, 3.06177191]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[3] == approx(s2[0])
+    assert y2[3] == approx(s2[1])
+
+
+def test_restart_celi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using CELI."""
+
+    op = dummy_operator.DummyOperator()
+    output_dir = "test_restart_celi"
+    op.output_dir = output_dir
+
+    # Perform simulation
+    dt = [0.75]
+    power = 1.0
+    openmc.deplete.celi(op, dt, power, print_out=False)
+
+    # Load the files
+    prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    # Re-create depletion operator and load previous results
+    op = dummy_operator.DummyOperator(prev_res)
+    op.output_dir = output_dir
+
+    # Perform restarts simulation
+    openmc.deplete.celi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [1.82078767, 0.97122898]
+    s2 = [2.68441779, 0.05125966]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[3] == approx(s2[0])
+    assert y2[3] == approx(s2[1])
+
+
+def test_restart_leqi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using LEQI."""
+
+    op = dummy_operator.DummyOperator()
+    output_dir = "test_restart_leqi"
+    op.output_dir = output_dir
+
+    # Perform simulation
+    dt = [0.75]
+    power = 1.0
+    openmc.deplete.leqi(op, dt, power, print_out=False)
+
+    # Load the files
+    prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    # Re-create depletion operator and load previous results
+    op = dummy_operator.DummyOperator(prev_res)
+    op.output_dir = output_dir
+
+    # Perform restarts simulation
+    openmc.deplete.leqi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [1.82078767, 0.97122898]
+    s2 = [2.74526197, 0.23339915]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[3] == approx(s2[0])
+    assert y2[3] == approx(s2[1])
+
+def test_restart_si_celi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using SI-CELI."""
+
+    op = dummy_operator.DummyOperator()
+    output_dir = "test_restart_si_celi"
+    op.output_dir = output_dir
+
+    # Perform simulation
+    dt = [0.75]
+    power = 1.0
+    openmc.deplete.si_celi(op, dt, power, print_out=False)
+
+    # Load the files
+    prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    # Re-create depletion operator and load previous results
+    op = dummy_operator.DummyOperator(prev_res)
+    op.output_dir = output_dir
+
+    # Perform restarts simulation
+    openmc.deplete.si_celi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.03325094, 1.16826254]
+    s2 = [2.69291933, 0.37907772]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[3] == approx(s2[0])
+    assert y2[3] == approx(s2[1])
+
+
+def test_restart_si_leqi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using SI-LEQI."""
+
+    op = dummy_operator.DummyOperator()
+    output_dir = "test_restart_si_leqi"
+    op.output_dir = output_dir
+
+    # Perform simulation
+    dt = [0.75]
+    power = 1.0
+    openmc.deplete.si_leqi(op, dt, power, print_out=False)
+
+    # Load the files
+    prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    # Re-create depletion operator and load previous results
+    op = dummy_operator.DummyOperator(prev_res)
+    op.output_dir = output_dir
+
+    # Perform restarts simulation
+    openmc.deplete.si_leqi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.03325094, 1.16826254]
+    s2 = [2.92711288, 0.53753236]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[3] == approx(s2[0])
+    assert y2[3] == approx(s2[1])
diff --git a/tests/unit_tests/test_deplete_si_celi.py b/tests/unit_tests/test_deplete_si_celi.py
new file mode 100644
index 000000000..8e84bcba8
--- /dev/null
+++ b/tests/unit_tests/test_deplete_si_celi.py
@@ -0,0 +1,37 @@
+"""Regression tests for openmc.deplete.integrator.si_celi algorithm.
+
+These tests integrate a simple test problem described in dummy_geometry.py.
+"""
+
+from pytest import approx
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+def test_si_celi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using si_celi"""
+
+    op = dummy_operator.DummyOperator()
+    op.output_dir = "test_integrator_regression"
+
+    # Perform simulation using the si_celi algorithm
+    dt = [0.75, 0.75]
+    power = 1.0
+    openmc.deplete.si_celi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.03325094, 1.16826254]
+    s2 = [2.69291933, 0.37907772]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[2] == approx(s2[0])
+    assert y2[2] == approx(s2[1])
diff --git a/tests/unit_tests/test_deplete_si_leqi.py b/tests/unit_tests/test_deplete_si_leqi.py
new file mode 100644
index 000000000..4814d8ab8
--- /dev/null
+++ b/tests/unit_tests/test_deplete_si_leqi.py
@@ -0,0 +1,37 @@
+"""Regression tests for openmc.deplete.integrator.si_leqi algorithm.
+
+These tests integrate a simple test problem described in dummy_geometry.py.
+"""
+
+from pytest import approx
+import openmc.deplete
+
+from tests import dummy_operator
+
+
+def test_si_leqi(run_in_tmpdir):
+    """Integral regression test of integrator algorithm using si_leqi"""
+
+    op = dummy_operator.DummyOperator()
+    op.output_dir = "test_integrator_regression"
+
+    # Perform simulation using the si_leqi algorithm
+    dt = [0.75, 0.75]
+    power = 1.0
+    openmc.deplete.si_leqi(op, dt, power, print_out=False)
+
+    # Load the files
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
+
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
+
+    # Reference solution
+    s1 = [2.03325094, 1.16826254]
+    s2 = [2.92711288, 0.53753236]
+
+    assert y1[1] == approx(s1[0])
+    assert y2[1] == approx(s1[1])
+
+    assert y1[2] == approx(s2[0])
+    assert y2[2] == approx(s2[1])