ConstraintSet::SubtractJacobianTerms allocates jac_buffer (and resizes concentrations for vectorized dense) on every call. This is both a performance hit inside the solver loop and can throw exceptions (bad_alloc), which is especially problematic since Rosenbrock::Solve is marked noexcept. Consider storing these scratch buffers as members of ConstraintSet or solver temporary variables and resizing once during setup.