diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py
index c23abef631..37269b7bd6 100644
--- a/tests/dummy_operator.py
+++ b/tests/dummy_operator.py
@@ -1,9 +1,76 @@
+from collections import namedtuple
+
 import numpy as np
 import scipy.sparse as sp
 from uncertainties import ufloat
 
 from openmc.deplete.reaction_rates import ReactionRates
 from openmc.deplete.abc import TransportOperator, OperatorResult
+from openmc.deplete import (
+    CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator,
+    EPCRK4Integrator, CF4Integrator, SICELIIntegrator, SILEQIIntegrator
+)
+
+# Bundle for nicely passing test data to depletion unit tests
+# solver should be a concrete subclass of openmc.deplete.abc.Integrator
+# atoms_1 should be the number of atoms of type 1 through the simulation
+# similar for atoms_2, but for type 2. This includes the first step
+# Solutions should be the exact solution that can be obtained using
+# the DummyOperator depletion matrix with two 0.75 second time steps
+DepletionSolutionTuple = namedtuple(
+    "DepletionSolutionTuple", "solver atoms_1 atoms_2")
+
+
+predictor_solution = DepletionSolutionTuple(
+    PredictorIntegrator, np.array([1.0, 2.46847546272295, 4.11525874568034]),
+    np.array([1.0, 0.986431226850467, -0.0581692232513460]))
+
+
+cecm_solution = DepletionSolutionTuple(
+    CECMIntegrator, np.array([1.0, 1.86872629872102, 2.18097439443550]),
+    np.array([1.0, 1.395525772416039, 2.69429754646747]))
+
+
+cf4_solution = DepletionSolutionTuple(
+    CF4Integrator, np.array([1.0, 2.06101629, 2.57241318]),
+    np.array([1.0, 1.37783588, 2.63731630]))
+
+
+epc_rk4_solution = DepletionSolutionTuple(
+    EPCRK4Integrator, np.array([1.0, 2.01978516, 2.05246421]),
+    np.array([1.0, 1.42038037, 3.06177191]))
+
+
+celi_solution = DepletionSolutionTuple(
+    CELIIntegrator, np.array([1.0, 1.82078767, 2.68441779]),
+    np.array([1.0, 0.97122898, 0.05125966]))
+
+
+si_celi_solution = DepletionSolutionTuple(
+    SICELIIntegrator, np.array([1.0, 2.03325094, 2.69291933]),
+    np.array([1.0, 1.16826254, 0.37907772]))
+
+
+leqi_solution = DepletionSolutionTuple(
+    LEQIIntegrator, np.array([1.0, 1.82078767, 2.74526197]),
+    np.array([1.0, 0.97122898, 0.23339915]))
+
+
+si_leqi_solution = DepletionSolutionTuple(
+    SILEQIIntegrator, np.array([1.0, 2.03325094, 2.92711288]),
+    np.array([1.0, 1.16826254, 0.53753236]))
+
+
+SCHEMES = {
+    "predictor": predictor_solution,
+    "cecm": cecm_solution,
+    "celi": celi_solution,
+    "cf4": cf4_solution,
+    "epc_rk4": epc_rk4_solution,
+    "leqi": leqi_solution,
+    "si_leqi": si_leqi_solution,
+    "si_celi": si_celi_solution,
+}
 
 
 class DummyOperator(TransportOperator):
@@ -21,6 +88,7 @@ class DummyOperator(TransportOperator):
     """
     def __init__(self, previous_results=None):
         self.prev_res = previous_results
+        self.output_dir = "."
 
     def __call__(self, vec, power, print_out=False):
         """Evaluates F(y)
@@ -76,7 +144,6 @@ class DummyOperator(TransportOperator):
         y_1 = rates[0, 0]
         y_2 = rates[1, 0]
 
-        mat = np.zeros((2, 2))
         a11 = np.sin(y_2)
         a12 = np.cos(y_1)
         a21 = -np.cos(y_2)
@@ -106,7 +173,8 @@ class DummyOperator(TransportOperator):
     def local_mats(self):
         """
         local_mats : list of str
-            A list of all material IDs to be burned.  Used for sorting the simulation.
+            A list of all material IDs to be burned.  Used for sorting the
+            simulation.
         """
 
         return ["1"]
@@ -153,7 +221,8 @@ class DummyOperator(TransportOperator):
         nuc_list : list of str
             A list of all nuclide names. Used for sorting the simulation.
         burn_list : list of int
-            A list of all cell IDs to be burned.  Used for sorting the simulation.
+            A list of all cell IDs to be burned.  Used for sorting the
+            simulation.
         full_burn_list : OrderedDict of str to int
             Maps cell name to index in global geometry.
 
diff --git a/tests/unit_tests/test_deplete_cecm.py b/tests/unit_tests/test_deplete_cecm.py
deleted file mode 100644
index 65dd61643b..0000000000
--- a/tests/unit_tests/test_deplete_cecm.py
+++ /dev/null
@@ -1,44 +0,0 @@
-"""Regression tests for openmc.deplete.integrator.cecm algorithm.
-
-These tests integrate a simple test problem described in dummy_geometry.py.
-"""
-
-from pytest import approx
-from openmc.deplete import CECMIntegrator, ResultsList
-
-from tests import dummy_operator
-
-
-def test_cecm(run_in_tmpdir):
-    """Integral regression test of integrator algorithm using CE/CM."""
-
-    op = dummy_operator.DummyOperator()
-    op.output_dir = "test_integrator_regression"
-
-    # Perform simulation using the MCNPX/MCNP6 algorithm
-    dt = [0.75, 0.75]
-    power = 1.0
-    integrator = CECMIntegrator(op, dt, power)
-    integrator.integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
-
-    _, y1 = res.get_atoms("1", "1")
-    _, y2 = res.get_atoms("1", "2")
-
-    # Mathematica solution
-    s1 = [1.86872629872102, 1.395525772416039]
-    s2 = [2.18097439443550, 2.69429754646747]
-
-    assert y1[1] == approx(s1[0])
-    assert y2[1] == approx(s1[1])
-
-    assert y1[2] == approx(s2[0])
-    assert y2[2] == approx(s2[1])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_celi.py b/tests/unit_tests/test_deplete_celi.py
deleted file mode 100644
index 01ec60d372..0000000000
--- a/tests/unit_tests/test_deplete_celi.py
+++ /dev/null
@@ -1,43 +0,0 @@
-"""Regression tests for openmc.deplete.integrator.celi algorithm.
-
-These tests integrate a simple test problem described in dummy_geometry.py.
-"""
-
-from pytest import approx
-from openmc.deplete import ResultsList, CELIIntegrator
-
-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
-    CELIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(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])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_cf4.py b/tests/unit_tests/test_deplete_cf4.py
deleted file mode 100644
index 817dfb5e17..0000000000
--- a/tests/unit_tests/test_deplete_cf4.py
+++ /dev/null
@@ -1,43 +0,0 @@
-"""Regression tests for openmc.deplete.integrator.cf4 algorithm.
-
-These tests integrate a simple test problem described in dummy_geometry.py.
-"""
-
-from pytest import approx
-from openmc.deplete import CF4Integrator, ResultsList
-
-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
-    CF4Integrator(op, dt, power).integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(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])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_epc_rk4.py b/tests/unit_tests/test_deplete_epc_rk4.py
deleted file mode 100644
index 0c656d7dbe..0000000000
--- a/tests/unit_tests/test_deplete_epc_rk4.py
+++ /dev/null
@@ -1,43 +0,0 @@
-"""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
-from openmc.deplete import EPCRK4Integrator, ResultsList
-
-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
-    EPCRK4Integrator(op, dt, power).integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(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])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index fd4e2ba3ff..47c9690ae9 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -7,7 +7,6 @@ will be left unimplemented and testing will be done via regression.
 """
 
 import copy
-import os
 from unittest.mock import MagicMock
 
 import numpy as np
@@ -18,6 +17,8 @@ from openmc.deplete import (
     ReactionRates, Results, ResultsList, comm, OperatorResult,
     PredictorIntegrator, SICELIIntegrator)
 
+from tests import dummy_operator
+
 
 def test_results_save(run_in_tmpdir):
     """Test data save module"""
@@ -41,10 +42,11 @@ def test_results_save(run_in_tmpdir):
         full_burn_list.append(str(2*i))
         full_burn_list.append(str(2*i + 1))
 
-    burn_list = full_burn_list[2*comm.rank : 2*comm.rank + 2]
+    burn_list = full_burn_list[2*comm.rank: 2*comm.rank + 2]
     nuc_list = ["na", "nb"]
 
-    op.get_results_info.return_value = vol_dict, nuc_list, burn_list, full_burn_list
+    op.get_results_info.return_value = (
+        vol_dict, nuc_list, burn_list, full_burn_list)
 
     # Construct x
     x1 = []
@@ -130,3 +132,30 @@ def test_bad_integrator_inputs(timesteps):
 
     with pytest.raises(ValueError, match="n_steps"):
         SICELIIntegrator(op, timesteps, [1], n_steps=0)
+
+
+@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
+def test_integrator(run_in_tmpdir, scheme):
+    """Test the integrators against their expected values"""
+
+    bundle = dummy_operator.SCHEMES[scheme]
+    operator = dummy_operator.DummyOperator()
+    bundle.solver(operator, [0.75, 0.75], 1.0).integrate()
+
+    # get expected results
+
+    res = ResultsList.from_hdf5(
+        operator.output_dir / "depletion_results.h5")
+
+    t1, y1 = res.get_atoms("1", "1")
+    t2, y2 = res.get_atoms("1", "2")
+
+    assert (t1 == [0.0, 0.75, 1.5]).all()
+    assert y1 == pytest.approx(bundle.atoms_1)
+    assert (t2 == [0.0, 0.75, 1.5]).all()
+    assert y2 == pytest.approx(bundle.atoms_2)
+
+    # test structure of depletion time dataset
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (2, )
+    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_leqi.py b/tests/unit_tests/test_deplete_leqi.py
deleted file mode 100644
index 113cb6e306..0000000000
--- a/tests/unit_tests/test_deplete_leqi.py
+++ /dev/null
@@ -1,43 +0,0 @@
-"""Regression tests for openmc.deplete.integrator.leqi algorithm.
-
-These tests integrate a simple test problem described in dummy_geometry.py.
-"""
-
-from pytest import approx
-from openmc.deplete import LEQIIntegrator, ResultsList
-
-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
-    LEQIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(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])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_predictor.py b/tests/unit_tests/test_deplete_predictor.py
deleted file mode 100644
index 38bee86c39..0000000000
--- a/tests/unit_tests/test_deplete_predictor.py
+++ /dev/null
@@ -1,44 +0,0 @@
-"""Regression tests for openmc.deplete.integrator.predictor algorithm.
-
-These tests integrate a simple test problem described in dummy_geometry.py.
-"""
-
-from pytest import approx
-from openmc.deplete import PredictorIntegrator, ResultsList
-
-from tests import dummy_operator
-
-
-def test_predictor(run_in_tmpdir):
-    """Integral regression test of integrator algorithm using predictor"""
-
-    op = dummy_operator.DummyOperator()
-    op.output_dir = "test_integrator_regression"
-
-    # Perform simulation using the predictor algorithm
-    dt = [0.75, 0.75]
-    power = 1.0
-    PredictorIntegrator(op, dt, power).integrate()
-
-
-    # Load the files
-    res = ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
-
-    _, y1 = res.get_atoms("1", "1")
-    _, y2 = res.get_atoms("1", "2")
-
-    # Mathematica solution
-    s1 = [2.46847546272295, 0.986431226850467]
-    s2 = [4.11525874568034, -0.0581692232513460]
-
-    assert y1[1] == approx(s1[0])
-    assert y2[1] == approx(s1[1])
-
-    assert y1[2] == approx(s2[0])
-    assert y2[2] == approx(s2[1])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py
index c55f7fd0f4..82ef05ac7b 100644
--- a/tests/unit_tests/test_deplete_restart.py
+++ b/tests/unit_tests/test_deplete_restart.py
@@ -4,96 +4,13 @@ These tests run in two steps, a first run then a restart run, a simple test
 problem described in dummy_geometry.py.
 """
 
-from pytest import approx, raises
+import pytest
+
 import openmc.deplete
-from openmc.deplete import (
-    CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator,
-    EPCRK4Integrator, CF4Integrator, SICELIIntegrator, SILEQIIntegrator
-)
 
 from tests import dummy_operator
 
 
-def test_restart_predictor(run_in_tmpdir):
-    """Integral regression test of integrator algorithm using predictor."""
-
-    op = dummy_operator.DummyOperator()
-    output_dir = "test_restart_predictor"
-    op.output_dir = output_dir
-
-    # Perform simulation using the predictor algorithm
-    dt = [0.75]
-    power = 1.0
-    PredictorIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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 using the predictor algorithm
-    PredictorIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
-
-    _, y1 = res.get_atoms("1", "1")
-    _, y2 = res.get_atoms("1", "2")
-
-    # Mathematica solution
-    s1 = [2.46847546272295, 0.986431226850467]
-    s2 = [4.11525874568034, -0.0581692232513460]
-
-    assert y1[1] == approx(s1[0])
-    assert y2[1] == approx(s1[1])
-
-    assert y1[2] == approx(s2[0])
-    assert y2[2] == approx(s2[1])
-
-
-def test_restart_cecm(run_in_tmpdir):
-    """Integral regression test of integrator algorithm using CE/CM."""
-
-    op = dummy_operator.DummyOperator()
-    output_dir = "test_restart_cecm"
-    op.output_dir = output_dir
-
-    # Perform simulation using the MCNPX/MCNP6 algorithm
-    dt = [0.75]
-    power = 1.0
-    cecm = CECMIntegrator(op, dt, power)
-    cecm.integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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 using the MCNPX/MCNP6 algorithm
-    cecm_restart = CECMIntegrator(op, dt, power)
-    cecm_restart.integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
-
-    _, y1 = res.get_atoms("1", "1")
-    _, y2 = res.get_atoms("1", "2")
-
-    # Mathematica solution
-    s1 = [1.86872629872102, 1.395525772416039]
-    s2 = [2.18097439443550, 2.69429754646747]
-
-    assert y1[1] == approx(s1[0])
-    assert y2[1] == approx(s1[1])
-
-    assert y1[2] == approx(s2[0])
-    assert y2[2] == approx(s2[1])
-
-
 def test_restart_predictor_cecm(run_in_tmpdir):
     """Test to ensure that schemes with different stages are not compatible"""
 
@@ -104,18 +21,19 @@ def test_restart_predictor_cecm(run_in_tmpdir):
     # Perform simulation using the predictor algorithm
     dt = [0.75]
     power = 1.0
-    PredictorIntegrator(op, dt, power).integrate()
+    openmc.deplete.PredictorIntegrator(op, dt, power).integrate()
 
     # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
+    prev_res = openmc.deplete.ResultsList.from_hdf5(
+        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
 
     # check ValueError is raised, indicating previous and current stages
-    with raises(ValueError, match="incompatible.* 1.*2"):
-        CECMIntegrator(op, dt, power)
+    with pytest.raises(ValueError, match="incompatible.* 1.*2"):
+        openmc.deplete.CECMIntegrator(op, dt, power)
 
 
 def test_restart_cecm_predictor(run_in_tmpdir):
@@ -129,249 +47,48 @@ def test_restart_cecm_predictor(run_in_tmpdir):
     # Perform simulation using the MCNPX/MCNP6 algorithm
     dt = [0.75]
     power = 1.0
-    cecm = CECMIntegrator(op, dt, power)
+    cecm = openmc.deplete.CECMIntegrator(op, dt, power)
     cecm.integrate()
 
     # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
+    prev_res = openmc.deplete.ResultsList.from_hdf5(
+        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
 
     # check ValueError is raised, indicating previous and current stages
-    with raises(ValueError, match="incompatible.* 2.*1"):
-        PredictorIntegrator(op, dt, power)
+    with pytest.raises(ValueError, match="incompatible.* 2.*1"):
+        openmc.deplete.PredictorIntegrator(op, dt, power)
 
-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
+@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
+def test_restart(run_in_tmpdir, scheme):
+    # set up the problem
 
-    # Perform simulation
-    dt = [0.75]
-    power = 1.0
-    CF4Integrator(op, dt, power).integrate()
+    bundle = dummy_operator.SCHEMES[scheme]
 
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
+    operator = dummy_operator.DummyOperator()
 
-    # Re-create depletion operator and load previous results
-    op = dummy_operator.DummyOperator(prev_res)
-    op.output_dir = output_dir
+    # take first step
+    bundle.solver(operator, [0.75], 1.0).integrate()
 
-    # Perform restarts simulation
-    CF4Integrator(op, dt, power).integrate()
+    # restart
+    prev_res = openmc.deplete.ResultsList.from_hdf5(
+        operator.output_dir / "depletion_results.h5")
+    operator = dummy_operator.DummyOperator(prev_res)
 
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
+    # take second step
+    bundle.solver(operator, [0.75], 1.0).integrate()
 
-    _, y1 = res.get_atoms("1", "1")
-    _, y2 = res.get_atoms("1", "2")
+    # compare results
 
-    # Reference solution
-    s1 = [2.06101629, 1.37783588]
-    s2 = [2.57241318, 2.63731630]
+    results = openmc.deplete.ResultsList.from_hdf5(
+        operator.output_dir / "depletion_results.h5")
 
-    assert y1[1] == approx(s1[0])
-    assert y2[1] == approx(s1[1])
+    _t, y1 = results.get_atoms("1", "1")
+    _t, y2 = results.get_atoms("1", "2")
 
-    assert y1[2] == approx(s2[0])
-    assert y2[2] == 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
-    EPCRK4Integrator(op, dt, power).integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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
-    EPCRK4Integrator(op, dt, power).integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(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])
-
-
-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
-    CELIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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
-    CELIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(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])
-
-
-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
-    LEQIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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
-    LEQIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(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])
-
-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
-    SICELIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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
-    SICELIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(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])
-
-
-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
-    nstages = 10
-    SILEQIIntegrator(op, dt, power, nstages).integrate()
-
-    # Load the files
-    prev_res = openmc.deplete.ResultsList.from_hdf5(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
-    SILEQIIntegrator(op, dt, power, nstages).integrate()
-
-    # Load the files
-    res = openmc.deplete.ResultsList.from_hdf5(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])
+    assert y1 == pytest.approx(bundle.atoms_1)
+    assert y2 == pytest.approx(bundle.atoms_2)
diff --git a/tests/unit_tests/test_deplete_si_celi.py b/tests/unit_tests/test_deplete_si_celi.py
deleted file mode 100644
index 5cc0d4e6c4..0000000000
--- a/tests/unit_tests/test_deplete_si_celi.py
+++ /dev/null
@@ -1,43 +0,0 @@
-"""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
-from openmc.deplete import SICELIIntegrator, ResultsList
-
-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
-    SICELIIntegrator(op, dt, power).integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(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])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)
diff --git a/tests/unit_tests/test_deplete_si_leqi.py b/tests/unit_tests/test_deplete_si_leqi.py
deleted file mode 100644
index 8e3c3a4648..0000000000
--- a/tests/unit_tests/test_deplete_si_leqi.py
+++ /dev/null
@@ -1,43 +0,0 @@
-"""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
-from openmc.deplete import SILEQIIntegrator, ResultsList
-
-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
-    SILEQIIntegrator(op, dt, power, 10).integrate()
-
-    # Load the files
-    res = ResultsList.from_hdf5(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])
-
-    # Test structure of depletion time dataset
-
-    dep_time = res.get_depletion_time()
-    assert dep_time.shape == (len(dt), )
-    assert all(dep_time > 0)