ODE solvers currently use std::vector<Polytope*> to represent domains, where a NULL pointer is used to indicate the unconstrained domain R^d.
This mixes the domain semantic with pointer validity and makes solver behavior hard to reason. In several places also solver steps proceed differently depending on whether this is a NULL pointer without exposing this information to the caller.
We should avoid using NULL as a sentinel, instead we should make the domain representation explicit. The handling of the unconstrained domain should live inside the corresponding polytope or domain abstraction not in the solver logic.