diff --git a/dwave/optimization/generators.py b/dwave/optimization/generators.py
index e95d9624..bfbb4d9a 100644
--- a/dwave/optimization/generators.py
+++ b/dwave/optimization/generators.py
@@ -34,6 +34,7 @@
     concatenate,
     exp,
     expit,
+    is_disjoint_cover,
     logical_or,
     maximum,
     minimum,
@@ -377,7 +378,7 @@ class as the decision variable being optimized, with permutations of its
         The :meth:`~dwave.optimization.model.Model.iter_decisions` method
         obtains the decision variables of the generated model.
 
-        >>> routes = next(model.iter_decisions())
+        >>> routes = list(model.iter_decisions())
 
         To test the solution above, set it in the model as the state of the
         decision variable. **Skip these next lines** if you have submitted your
@@ -385,7 +386,8 @@ class as the decision variable being optimized, with permutations of its
         nonlinear :term:`solver`.
 
         >>> model.states.resize(1)
-        >>> routes.set_state(0, [[2., 7., 1., 5.], [4., 3., 8., 6., 0.]])
+        >>> routes[0].set_state(0, [2., 7., 1., 5.])
+        >>> routes[1].set_state(0, [4., 3., 8., 6., 0.])
 
         You can use the :meth:`~dwave.optimization.model.Model.iter_constraints`
         method to check feasibility of constructed or returned solutions. Here,
@@ -398,7 +400,7 @@ class as the decision variable being optimized, with permutations of its
         ...     for i in range(model.states.size()):
         ...         if capacity_constraint.state(i):    # Filter on feasibility
         ...             print((f"Objective value #{i} is {model.objective.state(i).round(2)} for routes\n"
-        ...                    f"    {[r.state(i).tolist() for r in routes.iter_successors()]}"))
+        ...                    f"    {[r.state(i).tolist() for r in routes]}"))
         Objective value #0 is 423.8 for routes
             [[2.0, 7.0, 1.0, 5.0], [4.0, 3.0, 8.0, 6.0, 0.0]]
     """
@@ -511,10 +513,8 @@ class as the decision variable being optimized, with permutations of its
     demand = model.constant(customer_demand)
     capacity = model.constant(vehicle_capacity)
 
-    # Add the decision variable
-    routes = model.disjoint_lists_symbol(
-        primary_set_size=num_customers,
-        num_disjoint_lists=number_of_vehicles)
+    routes = [model.list(num_customers, min_size=0) for _ in range(number_of_vehicles)]
+    model.add_constraint(is_disjoint_cover(routes))
 
     # The objective is to minimize the distance traveled.
     # This is calculated by adding the distance from the depot to the 1st customer
@@ -631,7 +631,7 @@ class as the decision variable being optimized, with permutations of its
         The :meth:`~dwave.optimization.model.Model.iter_decisions` method
         obtains the decision variables of the generated model.
 
-        >>> routes = next(model.iter_decisions())
+        >>> routes = list(model.iter_decisions())
 
         To test the solution above, set it in the model as the state of the
         decision variable. **Skip these next lines** if you have submitted your
@@ -639,7 +639,8 @@ class as the decision variable being optimized, with permutations of its
         nonlinear :term:`solver`.
 
         >>> model.states.resize(1)
-        >>> routes.set_state(0, [[0, 2], [1]])
+        >>> routes[0].set_state(0, [0, 2])
+        >>> routes[1].set_state(0, [1])
 
         You can use the :meth:`~dwave.optimization.model.Model.iter_constraints`
         method to check feasibility of constructed or returned solutions. Here,
@@ -757,9 +758,8 @@ class as the decision variable being optimized, with permutations of its
     one = model.constant(1)
 
     # Add the decision variable
-    routes = model.disjoint_lists_symbol(
-        primary_set_size=num_customers,
-        num_disjoint_lists=number_of_vehicles)
+    routes = [model.list(num_customers, min_size=0) for _ in range(number_of_vehicles)]
+    model.add_constraint(is_disjoint_cover(routes))
 
     # Capacity constraint
     capacity_constraints = [(demand[routes[vehicle_idx]].sum() <= capacity)
diff --git a/dwave/optimization/include/dwave-optimization/nodes.hpp b/dwave/optimization/include/dwave-optimization/nodes.hpp
index 01efa8b5..7ebfa47b 100644
--- a/dwave/optimization/include/dwave-optimization/nodes.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes.hpp
@@ -28,5 +28,6 @@
 #include "dwave-optimization/nodes/numbers.hpp"
 #include "dwave-optimization/nodes/quadratic_model.hpp"
 #include "dwave-optimization/nodes/reduce.hpp"
+#include "dwave-optimization/nodes/set_routines.hpp"
 #include "dwave-optimization/nodes/testing.hpp"
 #include "dwave-optimization/nodes/unaryop.hpp"
diff --git a/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp b/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp
index 0ce7def6..e4272609 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp
@@ -22,6 +22,45 @@
 
 namespace dwave::optimization {
 
+class IsDisjointCoverNode : public ScalarOutputMixin<ArrayNode, false> {
+ public:
+    IsDisjointCoverNode(std::span<ArrayNode*> node_ptrs, ssize_t primary_set_size);
+
+    /// @copydoc Array::buff()
+    double const* buff(const State& state) const override;
+
+    /// @copydoc Node::commit()
+    void commit(State& state) const override;
+
+    /// @copydoc Array::diff()
+    std::span<const Update> diff(const State& state) const override;
+
+    /// @copydoc Node::initialize_state()
+    void initialize_state(State& state) const override;
+
+    /// @copydoc Array::integral()
+    bool integral() const override;
+
+    /// @copydoc Array::max()
+    double max() const override;
+
+    /// @copydoc Array::min()
+    double min() const override;
+
+    /// The size of the primary set to be covered
+    ssize_t primary_set_size() const { return primary_set_size_; };
+
+    /// @copydoc Node::propagate()
+    void propagate(State& state) const override;
+
+    /// @copydoc Node::revert()
+    void revert(State& state) const override;
+
+ private:
+    ssize_t primary_set_size_;
+    std::vector<Array*> operands_;
+};
+
 class IsInNode : public ArrayOutputMixin<ArrayNode> {
  public:
     IsInNode(ArrayNode* element_ptr, ArrayNode* test_elements_ptr);
diff --git a/dwave/optimization/libcpp/nodes/set_routines.pxd b/dwave/optimization/libcpp/nodes/set_routines.pxd
index 10bd7a1f..64d089da 100644
--- a/dwave/optimization/libcpp/nodes/set_routines.pxd
+++ b/dwave/optimization/libcpp/nodes/set_routines.pxd
@@ -16,5 +16,8 @@ from dwave.optimization.libcpp.graph cimport ArrayNode
 
 
 cdef extern from "dwave-optimization/nodes/set_routines.hpp" namespace "dwave::optimization" nogil:
+    cdef cppclass IsDisjointCoverNode(ArrayNode):
+        Py_ssize_t primary_set_size() const
+
     cdef cppclass IsInNode(ArrayNode):
         pass
diff --git a/dwave/optimization/mathematical.py b/dwave/optimization/mathematical.py
index f0aac7c1..7ead06df 100644
--- a/dwave/optimization/mathematical.py
+++ b/dwave/optimization/mathematical.py
@@ -38,6 +38,7 @@
     Exp,
     Expit,
     Extract,
+    IsDisjointCover,
     IsIn,
     LessEqual,
     LinearProgram,
@@ -95,6 +96,7 @@
     "expit",
     "extract",
     "hstack",
+    "is_disjoint_cover",
     "isin",
     "less_equal",
     "linprog",
@@ -1110,6 +1112,52 @@ def hstack(arrays: collections.abc.Sequence[ArraySymbol]) -> ArraySymbol:
         return concatenate(arrays, 1)
 
 
+def is_disjoint_cover(subsets: list[ArraySymbol], *, primary_set_size: int|None = None) -> IsDisjointCover:
+    """Return whether the symbols are disjoint, set-like, and cover a set of integers.
+
+    Determines whether a collection of array symbols is disjoint and the union equals a fixed set.
+
+    Args:
+        subsets: List of array symbols to test whether they are disjoint and cover the primary set
+        primary_set_size: Number of elements in the primary set: {0, 1, ..., primary_set_size - 1}.
+            Must be non-negative. If `None`, primary set size is inferred from the common, finite,
+            nonnegative maximum value of the arrays.
+
+    Returns:
+        A scalar boolean-valued array symbol indicating whether the subsets are a disjoint cover
+
+    Examples:
+        >>> from dwave.optimization.model import Model
+        >>> from dwave.optimization.mathematical import is_disjoint_cover
+        ...
+        >>> model = Model()
+        >>> model.states.resize(1)
+        >>> subsets = []
+        >>> subsets.append(model.constant([0, 1]))
+        >>> subsets.append(model.constant([2, 3, 4]))
+        >>> subsets.append(model.constant([]))
+        >>> is_disjoint = is_disjoint_cover(subsets, primary_set_size=5)
+        >>> with model.lock():
+        ...     print(is_disjoint.state(0))
+        1.0
+
+    See Also:
+        :class:`~dwave.optimization.symbols.IsDisjointCover`: Generated symbol
+
+        :func:`.isin`
+
+    .. versionadded:: 0.7.2
+    """
+    if primary_set_size is None:
+        primary_set_size = int(subsets[0].info().max) + 1
+        if not np.isfinite(primary_set_size) or primary_set_size <= 0:
+            raise ValueError("Cannot infer primary set size from subsets")
+        for subset in subsets[1:]:
+            if int(subset.info().max) != primary_set_size - 1:
+                raise ValueError("Cannot infer primary set size from subsets")
+    return IsDisjointCover(subsets, primary_set_size)
+
+
 def isin(element: ArraySymbol, test_elements: ArraySymbol) -> IsIn:
     """Return which values of one array symbol are in another.
 
diff --git a/dwave/optimization/model.py b/dwave/optimization/model.py
index 2cfa787b..7bdcc85b 100644
--- a/dwave/optimization/model.py
+++ b/dwave/optimization/model.py
@@ -591,6 +591,15 @@ def disjoint_lists_symbol(
 
             :meth:`.iter_decisions`, :meth:`.iter_successors`
         """
+        warnings.warn(
+            "The use of Model.disjoint_lists_symbol() is deprecated "
+            "since dwave.optimization 0.7.2. Use\n"
+            "from dwave.optimization.mathematical import is_disjoint_cover\n"
+            "lists = [model.list(primary_set_size, min_size=0) for _ in range(num_disjoint_lists)]\n"
+            "model.add_constraint(is_disjoint_cover(primary_set_size, lists))",
+            DeprecationWarning,
+        )
+
         from dwave.optimization.symbols import DisjointLists, DisjointList  # avoid circular import
         disjoint_lists = DisjointLists(self, primary_set_size, num_disjoint_lists)
 
diff --git a/dwave/optimization/src/nodes/set_routines.cpp b/dwave/optimization/src/nodes/set_routines.cpp
index a15301a4..c6ca5fd9 100644
--- a/dwave/optimization/src/nodes/set_routines.cpp
+++ b/dwave/optimization/src/nodes/set_routines.cpp
@@ -20,6 +20,182 @@
 
 namespace dwave::optimization {
 
+// IsDisjointCoverNode *******************************************************************
+
+struct IsDisjointCoverNodeData : public NodeStateData {
+ private:
+    // simple wrapper around ssize_t with default value 1, indicating no violations
+    struct Count_ {
+        ssize_t value = 1;
+
+        Count_& operator=(const ssize_t& rhs) {
+            this->value = rhs;
+            return *this;
+        }
+
+        Count_& operator+=(const ssize_t& rhs) {
+            this->value += rhs;
+            return *this;
+        }
+
+        Count_& operator-=(const ssize_t& rhs) {
+            this->value -= rhs;
+            return *this;
+        }
+
+        bool operator==(const ssize_t& rhs) { return this->value == rhs; }
+    };
+
+ public:
+    IsDisjointCoverNodeData(std::vector<ssize_t>& count) : is_disjoint_cover(0, 0.0, 0.0) {
+        // Record whether the count is not equal to 1 in the count_violations map
+        for (ssize_t i = 0, stop = count.size(); i < stop; ++i) {
+            if (count[i] != 1) {
+                count_violations[i] = count[i];
+            }
+        }
+        is_disjoint_cover.value = static_cast<double>(count_violations.size() == 0);
+        is_disjoint_cover.old = is_disjoint_cover.value;
+    };
+
+    void propagate() {
+        // Incorporate changed elements diff into the count_violations map
+        for (const auto& element : elements_decremented) {
+            count_violations[element] -= 1;
+        }
+        for (const auto& element : elements_incremented) {
+            count_violations[element] += 1;
+        }
+        // Take out any non-violations
+        for (auto it = count_violations.begin(); it != count_violations.end();) {
+            if (it->second == 1) {
+                it = count_violations.erase(it);
+            } else {
+                ++it;
+            }
+        }
+        is_disjoint_cover.value = static_cast<double>(count_violations.size() == 0);
+    };
+
+    void commit() {
+        elements_decremented.clear();
+        elements_incremented.clear();
+        is_disjoint_cover.old = is_disjoint_cover.value;
+    };
+
+    void revert() {
+        for (const auto& element : elements_decremented) {
+            // Reverse the decrement
+            count_violations[element] += 1;
+        }
+        for (const auto& element : elements_incremented) {
+            // Reverse the increment
+            count_violations[element] -= 1;
+        }
+        // We'll leave non-violation entries in count_violations - to be cleared in the next
+        // propagate
+        elements_decremented.clear();
+        elements_incremented.clear();
+        is_disjoint_cover.value = is_disjoint_cover.old;
+    };
+
+    // The actual logical output; is_disjoint_cover.value is whether the current state is a disjoint
+    // cover
+    Update is_disjoint_cover;
+    // count_violations[i] = the number of times that element `i` appears in the predecessors, if
+    // not equal to 1
+    std::unordered_map<ssize_t, Count_> count_violations;
+    // The "diffs" - lists of elements that were removed (respectively, inserted) in the predecessor
+    // diffs
+    std::vector<ssize_t> elements_decremented;
+    std::vector<ssize_t> elements_incremented;
+};
+
+IsDisjointCoverNode::IsDisjointCoverNode(
+    std::span<ArrayNode*> node_ptrs,
+    ssize_t primary_set_size
+) :
+    ScalarOutputMixin<ArrayNode, false>(), primary_set_size_(primary_set_size) {
+    for (const auto& node : node_ptrs) {
+        if (!node->integral()) {
+            throw std::invalid_argument("Predecessors of DisjointCoverNode must be integral");
+        }
+        if (!(node->max() < primary_set_size)) {
+            throw std::invalid_argument("Predecessor exceeds primary set size");
+        }
+        if (!(node->min() >= 0)) {
+            throw std::invalid_argument("Predecessor exceeds primary set size");
+        }
+        add_predecessor_(node);
+        operands_.push_back(node);
+    }
+}
+
+void IsDisjointCoverNode::initialize_state(State& state) const {
+    std::vector<ssize_t> count(primary_set_size_, 0);
+    for (const auto& node : operands_) {
+        for (const auto& value : node->view(state)) {
+            ssize_t element = static_cast<ssize_t>(value);
+            // Should not be possible because of checks in constructor
+            assert(element < primary_set_size_);
+            assert(element >= 0);
+            count[element] += 1;
+        }
+    }
+    emplace_data_ptr_<IsDisjointCoverNodeData>(state, count);
+}
+
+void IsDisjointCoverNode::commit(State& state) const {
+    data_ptr_<IsDisjointCoverNodeData>(state)->commit();
+}
+
+void IsDisjointCoverNode::propagate(State& state) const {
+    IsDisjointCoverNodeData* data = data_ptr_<IsDisjointCoverNodeData>(state);
+
+    bool no_update = true;
+    for (const auto& pred : operands_) {
+        for (const auto& update : pred->diff(state)) {
+            no_update = false;
+            if (!update.placed()) {  // i.e. removed or changed.
+                auto element = static_cast<ssize_t>(update.old);
+                data->elements_decremented.push_back(element);
+            }
+            if (!update.removed()) {  // i.e. placed or changed.
+                auto element = static_cast<ssize_t>(update.value);
+                data->elements_incremented.push_back(element);
+            }
+        }
+    }
+
+    // If no update, this node is unchanged
+    if (no_update) {
+        return;
+    }
+    // Else, propagate the populated diffs to rest of the state
+    data->propagate();
+}
+
+void IsDisjointCoverNode::revert(State& state) const {
+    data_ptr_<IsDisjointCoverNodeData>(state)->revert();
+}
+
+double const* IsDisjointCoverNode::buff(const State& state) const {
+    return &data_ptr_<IsDisjointCoverNodeData>(state)->is_disjoint_cover.value;
+}
+
+std::span<const Update> IsDisjointCoverNode::diff(const State& state) const {
+    auto data = data_ptr_<IsDisjointCoverNodeData>(state);
+    return std::span(
+        &data->is_disjoint_cover, data->is_disjoint_cover.old != data->is_disjoint_cover.value
+    );
+}
+
+bool IsDisjointCoverNode::integral() const { return true; }
+
+double IsDisjointCoverNode::min() const { return 0.0; }
+
+double IsDisjointCoverNode::max() const { return 1.0; }
+
 // IsInNode *******************************************************************
 struct IsInNodeSetData {
     IsInNodeSetData() = default;
diff --git a/dwave/optimization/symbols/__init__.py b/dwave/optimization/symbols/__init__.py
index 478bf711..3ac14617 100644
--- a/dwave/optimization/symbols/__init__.py
+++ b/dwave/optimization/symbols/__init__.py
@@ -91,7 +91,7 @@
     Prod,
     Sum,
 )
-from dwave.optimization.symbols.set_routines import IsIn
+from dwave.optimization.symbols.set_routines import IsDisjointCover, IsIn
 from dwave.optimization.symbols.softmax import SoftMax
 from dwave.optimization.symbols.sorting import ArgSort
 from dwave.optimization.symbols.statistics import Mean
@@ -142,6 +142,7 @@
     "Extract",
     "Input",
     "IntegerVariable",
+    "IsDisjointCover",
     "IsIn",
     "LessEqual",
     "LinearProgram",
diff --git a/dwave/optimization/symbols/set_routines.pyi b/dwave/optimization/symbols/set_routines.pyi
index 83ff3422..9fed4c72 100644
--- a/dwave/optimization/symbols/set_routines.pyi
+++ b/dwave/optimization/symbols/set_routines.pyi
@@ -14,4 +14,6 @@
 
 from dwave.optimization.model import ArraySymbol as _ArraySymbol
 
+class IsDisjointCover(_ArraySymbol): ...
+
 class IsIn(_ArraySymbol): ...
diff --git a/dwave/optimization/symbols/set_routines.pyx b/dwave/optimization/symbols/set_routines.pyx
index bb57618e..3342969e 100644
--- a/dwave/optimization/symbols/set_routines.pyx
+++ b/dwave/optimization/symbols/set_routines.pyx
@@ -14,10 +14,75 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
+import collections.abc
+import json
+
 from cython.operator cimport typeid
+from libcpp.vector cimport vector
+
+
+from dwave.optimization._model cimport _Graph, _register, ArraySymbol, Symbol
+from dwave.optimization.libcpp cimport dynamic_cast_ptr
+from dwave.optimization.libcpp.graph cimport ArrayNode
+from dwave.optimization.libcpp.nodes.set_routines cimport IsDisjointCoverNode, IsInNode
+
+
+cdef class IsDisjointCover(ArraySymbol):
+    """Tests whether the symbols are disjoint, set-like, and cover a set of integers
+
+    See Also:
+
+    .. versionadded:: 0.7.2
+    """
+    def __init__(self, object inputs, Py_ssize_t n):
+        if (not isinstance(inputs, collections.abc.Sequence) or
+                not all(isinstance(arr, ArraySymbol) for arr in inputs)):
+            raise TypeError("disjoint_cover takes a sequence of array symbols")
+
+        if len(inputs) < 1:
+            raise ValueError("need at least one array symbol to to form a disjoint cover")
+
+        cdef _Graph model = inputs[0].model
+        cdef vector[ArrayNode*] cppinputs
+
+        for symbol in inputs:
+            if symbol.model is not model:
+                raise ValueError("all predecessors must be from the same model")
+            cppinputs.push_back((<ArraySymbol?>symbol).array_ptr)
+
+        self.primary_set_size = n
+        cdef IsDisjointCoverNode* ptr = model._graph.emplace_node[IsDisjointCoverNode](cppinputs, n)
+        self.initialize_arraynode(model, ptr)
+
+    @classmethod
+    def _from_symbol(cls, Symbol symbol):
+        cdef IsDisjointCoverNode* ptr = dynamic_cast_ptr[IsDisjointCoverNode](symbol.node_ptr)
+        if not ptr:
+            raise TypeError(f"given symbol cannot construct a {cls.__name__}")
+        cdef IsDisjointCover sym = cls.__new__(cls)
+        sym.primary_set_size = ptr.primary_set_size()
+        sym.initialize_arraynode(symbol.model, ptr)
+        return sym
+
+    @classmethod
+    def _from_zipfile(cls, zf, directory, _Graph model, predecessors):
+        with zf.open(directory + "args.json", "r") as f:
+            args = json.load(f)
+        return cls(list(predecessors), args["primary_set_size"])
+
+    def _into_zipfile(self, zf, directory):
+        super()._into_zipfile(zf, directory)
+
+        encoder = json.JSONEncoder(separators=(',', ':'))
+
+        # get the non-array args
+        args = dict()
+        args.update(primary_set_size=int(self.primary_set_size))
+        zf.writestr(directory + "args.json", encoder.encode(args))
+
+    cdef Py_ssize_t primary_set_size
 
-from dwave.optimization._model cimport _Graph, _register, ArraySymbol
-from dwave.optimization.libcpp.nodes.set_routines cimport IsInNode
+_register(IsDisjointCover, typeid(IsDisjointCoverNode))
 
 
 cdef class IsIn(ArraySymbol):
diff --git a/releasenotes/notes/is-disjoint-cover-9f6bd97b637511fc.yaml b/releasenotes/notes/is-disjoint-cover-9f6bd97b637511fc.yaml
new file mode 100644
index 00000000..c6615690
--- /dev/null
+++ b/releasenotes/notes/is-disjoint-cover-9f6bd97b637511fc.yaml
@@ -0,0 +1,11 @@
+---
+features:
+  - |
+    Add C++ ``IsDisjointCoverNode`` node/symbol, which computes whether
+    a collection of lists are disjoint and cover some primary set.
+  - Add ``IsDisjointCover`` symbol and `is_disjoint_cover` function.
+
+deprecations:
+  - |
+    Deprecates `disjoint_list_symbol`
+
diff --git a/tests/cpp/nodes/test_set_routines.cpp b/tests/cpp/nodes/test_set_routines.cpp
index 6770605f..93c1bb7b 100644
--- a/tests/cpp/nodes/test_set_routines.cpp
+++ b/tests/cpp/nodes/test_set_routines.cpp
@@ -28,6 +28,230 @@ using Catch::Matchers::RangeEquals;
 
 namespace dwave::optimization {
 
+TEST_CASE("IsDisjointCoverNode") {
+    auto graph = Graph();
+
+    GIVEN("Two constant nodes that are disjoint") {
+        auto c1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2});
+        auto c2_ptr = graph.emplace_node<ConstantNode>(std::vector{3, 4});
+        std::vector<ArrayNode*> sets{c1_ptr, c2_ptr};
+
+        THEN("We can construct a IsDisjointCover node") {
+            auto dc_ptr = graph.emplace_node<IsDisjointCoverNode>(sets, 5);
+            graph.emplace_node<ArrayValidationNode>(dc_ptr);
+
+            CHECK(dc_ptr->min() == 0.0);
+            CHECK(dc_ptr->max() == 1.0);
+            CHECK(dc_ptr->integral());
+            CHECK(dc_ptr->logical());
+
+            AND_WHEN("We initialize a state") {
+                auto state = graph.initialize_state();
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                }
+            }
+        }
+        THEN("We can construct a IsDisjointCover node with a larger primary set size") {
+            auto dc_ptr = graph.emplace_node<IsDisjointCoverNode>(sets, 6);
+            graph.emplace_node<ArrayValidationNode>(dc_ptr);
+
+            CHECK(dc_ptr->min() == 0.0);
+            CHECK(dc_ptr->max() == 1.0);
+            CHECK(dc_ptr->integral());
+            CHECK(dc_ptr->logical());
+
+            AND_WHEN("We initialize a state") {
+                auto state = graph.initialize_state();
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                }
+            }
+        }
+        THEN("We cannot construct a IsDisjointCover node with a smaller primary set size") {
+            CHECK_THROWS(graph.emplace_node<IsDisjointCoverNode>(sets, 4));
+        }
+    }
+
+    GIVEN("Two constant nodes that are disjoint but not sets") {
+        auto c1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 2});
+        auto c2_ptr = graph.emplace_node<ConstantNode>(std::vector{3, 4});
+        std::vector<ArrayNode*> sets{c1_ptr, c2_ptr};
+
+        THEN("We can construct a IsDisjointCover node") {
+            auto dc_ptr = graph.emplace_node<IsDisjointCoverNode>(sets, 5);
+            graph.emplace_node<ArrayValidationNode>(dc_ptr);
+
+            CHECK(dc_ptr->min() == 0.0);
+            CHECK(dc_ptr->max() == 1.0);
+            CHECK(dc_ptr->integral());
+            CHECK(dc_ptr->logical());
+
+            AND_WHEN("We initialize a state") {
+                auto state = graph.initialize_state();
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                }
+            }
+        }
+    }
+
+    GIVEN("Two set nodes") {
+        auto set1_ptr = graph.emplace_node<SetNode>(4);
+        auto set2_ptr = graph.emplace_node<SetNode>(5);
+        std::vector<ArrayNode*> sets{set1_ptr, set2_ptr};
+
+        THEN("We can construct a IsDisjointCover node") {
+            auto dc_ptr = graph.emplace_node<IsDisjointCoverNode>(sets, 5);
+            graph.emplace_node<ArrayValidationNode>(dc_ptr);
+
+            CHECK(dc_ptr->min() == 0.0);
+            CHECK(dc_ptr->max() == 1.0);
+            CHECK(dc_ptr->integral());
+            CHECK(dc_ptr->logical());
+
+            AND_WHEN("We initialize a state with disjoint sets that cover the primary set") {
+                auto state = graph.initialize_state();
+                set1_ptr->assign(state, std::vector<double>{0, 1, 2});
+                set2_ptr->assign(state, std::vector<double>{3, 4});
+                graph.propagate(state);
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                }
+
+                AND_WHEN("We commit and change the state") {
+                    graph.commit(state);
+                    set1_ptr->assign(state, std::vector<double>{0, 1, 2});
+                    set2_ptr->assign(state, std::vector<double>{3});
+                    graph.propagate(state);
+
+                    THEN("The IsDisjointCover state is correct") {
+                        CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                    }
+
+                    AND_WHEN("We revert") {
+                        graph.revert(state);
+
+                        THEN("The IsDisjointCover state is correct") {
+                            CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                        }
+                    }
+                }
+            }
+            AND_WHEN("We initialize a state with not-disjoint sets") {
+                auto state = graph.initialize_state();
+                set1_ptr->assign(state, std::vector<double>{0, 1, 2});
+                set2_ptr->assign(state, std::vector<double>{2, 3, 4});
+                graph.propagate(state);
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                }
+
+                AND_WHEN("We commit and change the state") {
+                    graph.commit(state);
+                    set1_ptr->assign(state, std::vector<double>{3, 1, 2});
+                    set2_ptr->assign(state, std::vector<double>{4, 0});
+                    graph.propagate(state);
+
+                    THEN("The IsDisjointCover state is correct") {
+                        CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                    }
+
+                    AND_WHEN("We revert") {
+                        graph.revert(state);
+
+                        THEN("The IsDisjointCover state is correct") {
+                            CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                        }
+                    }
+                }
+            }
+            AND_WHEN("We initialize a state with disjoint sets that do not cover the primary set") {
+                auto state = graph.initialize_state();
+                set1_ptr->assign(state, std::vector<double>{0, 1, 2});
+                set2_ptr->assign(state, std::vector<double>{3});
+                graph.propagate(state);
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                }
+
+                AND_WHEN("We commit and change the state") {
+                    graph.commit(state);
+                    set1_ptr->assign(state, std::vector<double>{0, 1, 3});
+                    set2_ptr->assign(state, std::vector<double>{2, 4});
+                    graph.propagate(state);
+
+                    THEN("The IsDisjointCover state is correct") {
+                        CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                    }
+
+                    AND_WHEN("We revert") {
+                        graph.revert(state);
+
+                        THEN("The IsDisjointCover state is correct") {
+                            CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    GIVEN("Three list nodes") {
+        auto list1_ptr = graph.emplace_node<ListNode>(5, 0, 5);
+        auto list2_ptr = graph.emplace_node<ListNode>(5, 0, 5);
+        auto list3_ptr = graph.emplace_node<ListNode>(5, 0, 5);
+        std::vector<ArrayNode*> lists{list1_ptr, list2_ptr, list3_ptr};
+
+        THEN("We can construct a IsDisjointCover node") {
+            auto dc_ptr = graph.emplace_node<IsDisjointCoverNode>(lists, 5);
+            graph.emplace_node<ArrayValidationNode>(dc_ptr);
+
+            CHECK(dc_ptr->min() == 0.0);
+            CHECK(dc_ptr->max() == 1.0);
+            CHECK(dc_ptr->integral());
+            CHECK(dc_ptr->logical());
+
+            AND_WHEN("We initialize a state with disjoint lists that cover the primary set") {
+                auto state = graph.initialize_state();
+                list1_ptr->assign(state, std::vector<double>{0, 1, 2});
+                list2_ptr->assign(state, std::vector<double>{3, 4});
+                list3_ptr->assign(state, std::vector<double>{});
+                graph.propagate(state);
+
+                THEN("The initial IsDisjointCover state is correct") {
+                    CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                }
+
+                AND_WHEN("We commit and change the state") {
+                    graph.commit(state);
+                    list1_ptr->exchange(state, 0, 2);
+                    graph.propagate(state);
+
+                    THEN("The IsDisjointCover state is correct") {
+                        CHECK_THAT(dc_ptr->view(state), RangeEquals({1.0}));
+                    }
+                }
+                AND_WHEN("We commit and change the state") {
+                    graph.commit(state);
+                    list1_ptr->shrink(state);
+                    graph.propagate(state);
+
+                    THEN("The IsDisjointCover state is correct") {
+                        CHECK_THAT(dc_ptr->view(state), RangeEquals({0.0}));
+                    }
+                }
+            }
+        }
+    }
+}
+
 TEST_CASE("IsInNode") {
     auto graph = Graph();
 
@@ -231,39 +455,40 @@ TEST_CASE("IsInNode") {
     }
 
     GIVEN("Two dynamic set nodes and an isin node") {
-            auto set1_ptr = graph.emplace_node<SetNode>(6);
-            auto set2_ptr = graph.emplace_node<SetNode>(6);
-            auto isin_ptr = graph.emplace_node<IsInNode>(set1_ptr, set2_ptr);
-            graph.emplace_node<ArrayValidationNode>(isin_ptr);
+        auto set1_ptr = graph.emplace_node<SetNode>(6);
+        auto set2_ptr = graph.emplace_node<SetNode>(6);
+        auto isin_ptr = graph.emplace_node<IsInNode>(set1_ptr, set2_ptr);
+        graph.emplace_node<ArrayValidationNode>(isin_ptr);
 
-            CHECK(isin_ptr->size() == -1);
+        CHECK(isin_ptr->size() == -1);
 
-            WHEN("We initialize a state") {
-                auto state = graph.initialize_state();
+        WHEN("We initialize a state") {
+            auto state = graph.initialize_state();
 
-                AND_WHEN("We assign the set nodes and propagate") {
-                    set1_ptr->assign(state, std::vector<double>{2});
-                    set2_ptr->assign(state, std::vector<double>{2});
-                    graph.propagate(state);
+            AND_WHEN("We assign the set nodes and propagate") {
+                set1_ptr->assign(state, std::vector<double>{2});
+                set2_ptr->assign(state, std::vector<double>{2});
+                graph.propagate(state);
 
-                    THEN("The isin's state is correct") {
-                        CHECK_THAT(isin_ptr->view(state), RangeEquals({1.0}));
-                    }
+                THEN("The isin's state is correct") {
+                    CHECK_THAT(isin_ptr->view(state), RangeEquals({1.0}));
+                }
 
-                    AND_WHEN("We commit, shrink the state of both set nodes, and propagate") {
-                        graph.commit(state);
+                AND_WHEN("We commit, shrink the state of both set nodes, and propagate") {
+                    graph.commit(state);
 
-                        set1_ptr->shrink(state);
-                        set2_ptr->shrink(state);
+                    set1_ptr->shrink(state);
+                    set2_ptr->shrink(state);
 
-                        graph.propagate(state);
+                    graph.propagate(state);
 
-                        THEN("The isin's state is correct") {
-                            CHECK(isin_ptr->view(state).size() == 0);
-                        }
+                    THEN("The isin's state is correct") {
+                        CHECK(isin_ptr->view(state).size() == 0);
                     }
                 }
             }
         }
     }
+}
+
 }  // namespace dwave::optimization
diff --git a/tests/test_generators.py b/tests/test_generators.py
index b24864ad..7d37f2e3 100644
--- a/tests/test_generators.py
+++ b/tests/test_generators.py
@@ -269,16 +269,17 @@ def test_basics(self):
             locations_y=[2, 0],
             depot_x_y=[0, 0])
 
-        self.assertEqual(model.num_decisions(), 1)
-        self.assertEqual(model.num_constraints(), num_vehicles)
+        self.assertEqual(model.num_decisions(), 2)
+        self.assertEqual(model.num_constraints(), num_vehicles+1)
         self.assertEqual(model.num_nodes(), 34)
         self.assertEqual(model.num_edges(), 46)
         self.assertEqual(model.is_locked(), True)
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        self.assertEqual(route.num_disjoint_lists(), 2)
-        route.set_state(0, [[0], [1]])
+        routes = list(model.iter_decisions())
+        self.assertEqual(len(routes), 2)
+        routes[0].set_state(0, [0])
+        routes[1].set_state(0, [1])
         self.assertEqual(model.objective.state(0), 10)
 
         # Depot at (0, 1.5) from demand[0] == 0
@@ -290,9 +291,10 @@ def test_basics(self):
             locations_y=[1.5, 1.5])
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        self.assertEqual(route.num_disjoint_lists(), 2)
-        route.set_state(0, [[0], []])
+        routes = list(model.iter_decisions())
+        self.assertEqual(len(routes), 2)
+        routes[0].set_state(0, [0])
+        routes[1].set_state(0, [])
         self.assertEqual(model.objective.state(0), 6)
 
         # Depot at (3, 2) from depot_x_y
@@ -305,9 +307,10 @@ def test_basics(self):
             depot_x_y=[3, 2])
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        self.assertEqual(route.num_disjoint_lists(), 2)
-        route.set_state(0, [[0], [1]])
+        routes = list(model.iter_decisions())
+        self.assertEqual(len(routes), 2)
+        routes[0].set_state(0, [0])
+        routes[1].set_state(0, [1])
         self.assertEqual(model.objective.state(0), 10)
 
         # Test asymmetric distances
@@ -318,9 +321,10 @@ def test_basics(self):
             distances=[[0, 1, 3], [2, 0, np.sqrt(17)], [4, np.sqrt(13), 0]])
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        self.assertEqual(route.num_disjoint_lists(), 2)
-        route.set_state(0, [[0], [1]])
+        routes = list(model.iter_decisions())
+        self.assertEqual(len(routes), 2)
+        routes[0].set_state(0, [0])
+        routes[1].set_state(0, [1])
         self.assertEqual(model.objective.state(0), 10)
 
     def test_cvrplib_P_n19_k2(self):
@@ -333,9 +337,9 @@ def test_cvrplib_P_n19_k2(self):
             vehicle_capacity=160)
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        route.set_state(0, [[i - 1 for i in [4, 11, 14, 12, 3, 17, 16, 8, 6]],
-                            [i - 1 for i in [18, 5, 13, 15, 9, 7, 2, 10, 1]]])
+        routes = list(model.iter_decisions())
+        routes[0].set_state(0, [i - 1 for i in [4, 11, 14, 12, 3, 17, 16, 8, 6]])
+        routes[1].set_state(0, [i - 1 for i in [18, 5, 13, 15, 9, 7, 2, 10, 1]])
         self.assertGreater(model.objective.state(0), 212)
         self.assertLess(model.objective.state(0), 213)
 
@@ -366,10 +370,10 @@ def test_state_serialization(self):
 
         model.states.resize(1)
 
-        routes, = model.iter_decisions()
+        routes = list(model.iter_decisions())
 
-        routes.set_state(0, [[i - 1 for i in [4, 11, 14, 12, 3, 17, 16, 8, 6]],
-                             [i - 1 for i in [18, 5, 13, 15, 9, 7, 2, 10, 1]]])
+        routes[0].set_state(0, [i - 1 for i in [4, 11, 14, 12, 3, 17, 16, 8, 6]])
+        routes[1].set_state(0, [i - 1 for i in [18, 5, 13, 15, 9, 7, 2, 10, 1]])
 
         # just smoke test
         with model.states.to_file() as f:
@@ -526,14 +530,15 @@ def test_basics(self):
             service_time=[0, 0, 0])
 
         min_expected_number_of_constraints = num_vehicles * 2 + n_time_windows + n_customers
-        self.assertEqual(model.num_decisions(), 1)
+        self.assertEqual(model.num_decisions(), 2)
         self.assertGreaterEqual(model.num_constraints(),min_expected_number_of_constraints)
         self.assertEqual(model.is_locked(), True)
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        self.assertEqual(route.num_disjoint_lists(), 2)
-        route.set_state(0, [[0], [1]])
+        routes = list(model.iter_decisions())
+        self.assertEqual(len(routes), 2)
+        routes[0].set_state(0, [0])
+        routes[1].set_state(0, [1])
         self.assertEqual(model.objective.state(0), 15)
 
         # Test asymmetric distances
@@ -548,9 +553,10 @@ def test_basics(self):
         )
 
         model.states.resize(1)
-        route = next(model.iter_decisions())
-        self.assertEqual(route.num_disjoint_lists(), 2)
-        route.set_state(0, [[0], [1]])
+        routes = list(model.iter_decisions())
+        self.assertEqual(len(routes), 2)
+        routes[0].set_state(0, [0])
+        routes[1].set_state(0, [1])
         self.assertEqual(model.objective.state(0), 10)
 
     def test_serialization(self):
@@ -600,10 +606,10 @@ def test_state_serialization(self):
 
         model.states.resize(1)
 
-        routes, = model.iter_decisions()
+        routes = list(model.iter_decisions())
 
-        routes.set_state(0, [[i-1 for i in [1, 2]],
-                             [i-1 for i in [3]]])
+        routes[0].set_state(0, [i-1 for i in [1, 2]])
+        routes[1].set_state(0, [i-1 for i in [3]])
 
         # just smoke test
         with model.states.to_file() as f:
diff --git a/tests/test_model.py b/tests/test_model.py
index 41de0f34..62f898bf 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -379,7 +379,9 @@ def test_remove_unused_symbols(self):
             self.assertEqual(num_removed, 2)
             self.assertEqual(model.num_symbols(), 1)
 
+
         with self.subTest("disjoint lists"):
+            self.skipTest("Deprecated symbol")
             model = Model()
 
             disjoint_lists = model.disjoint_lists_symbol(10, 4)
diff --git a/tests/test_symbols.py b/tests/test_symbols.py
index 8c5bf18b..7ec57633 100644
--- a/tests/test_symbols.py
+++ b/tests/test_symbols.py
@@ -1497,6 +1497,52 @@ def test_state_serialization_explicit(self):
         np.testing.assert_array_equal(ys[2].state(), [0, 0, 0, 0, 1])
 
 
+class TestDisjointCover(utils.SymbolTests):
+    def generate_symbols(self):
+        model = Model()
+        sets = [
+            model.constant([0, 1, 2]),
+            model.constant([3, 4])
+        ]
+        cover = dwave.optimization.symbols.IsDisjointCover(sets, 5)
+
+        with model.lock():
+            yield cover
+
+    def test(self):
+        from dwave.optimization.symbols import IsDisjointCover
+        model = Model()
+        sets = [
+            model.constant([0, 1, 2]),
+            model.constant([3, 4])
+        ]
+        cover = dwave.optimization.symbols.IsDisjointCover(sets, 5)
+        self.assertIsInstance(cover, IsDisjointCover)
+
+    def test_state(self):
+        model = Model()
+        # a disjoint cover
+        sets = [
+            model.constant([0, 1, 2]),
+            model.constant([3, 4])
+        ]
+        cover = dwave.optimization.symbols.IsDisjointCover(sets, 5)
+        model.states.resize(1)
+        with model.lock():
+            expected = np.array([1.0])
+            np.testing.assert_array_almost_equal(cover.state(0), expected)
+
+        # not disjoint
+        sets = [
+            model.constant([0, 1, 2]),
+            model.constant([2, 3, 4])
+        ]
+        cover = dwave.optimization.symbols.IsDisjointCover(sets, 5)
+        with model.lock():
+            expected = np.array([0.0])
+            np.testing.assert_array_almost_equal(cover.state(0), expected)
+
+@unittest.skip("Deprecated symbol")
 class TestDisjointListsVariable(utils.SymbolTests):
     def test_inequality(self):
         # TODO re-enable this once equality has been fixed