diff --git a/.gitignore b/.gitignore
index 750d2605..5ec114b2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -79,6 +79,7 @@ subprojects/fourdst/
 subprojects/libplugin/
 subprojects/minizip-ng-4.0.10/
 subprojects/cvode-*/
+subprojects/kinsol-*/
 *.fbundle
 *.wraplock
 
diff --git a/build-config/meson.build b/build-config/meson.build
index 619b242e..7008fe82 100644
--- a/build-config/meson.build
+++ b/build-config/meson.build
@@ -3,7 +3,7 @@ cmake = import('cmake')
 subdir('fourdst')
 subdir('libplugin')
 
-subdir('cvode')
+subdir('sundials')
 
 subdir('boost')
 subdir('cppad')
diff --git a/build-config/cvode/meson.build b/build-config/sundials/cvode/meson.build
similarity index 97%
rename from build-config/cvode/meson.build
rename to build-config/sundials/cvode/meson.build
index c0a9e74d..cc382cb8 100644
--- a/build-config/cvode/meson.build
+++ b/build-config/sundials/cvode/meson.build
@@ -34,7 +34,7 @@ sundials_sunmatrixdense_dep = cvode_sp.dependency('sundials_sunmatrixdense_share
 # For the dense linear solver library
 sundials_sunlinsoldense_dep = cvode_sp.dependency('sundials_sunlinsoldense_shared')
 
-sundials_dep = declare_dependency(
+cvode_dep = declare_dependency(
     dependencies: [
         sundials_core_dep,
         sundials_cvode_dep,
diff --git a/build-config/sundials/kinsol/meson.build b/build-config/sundials/kinsol/meson.build
new file mode 100644
index 00000000..3c780841
--- /dev/null
+++ b/build-config/sundials/kinsol/meson.build
@@ -0,0 +1,29 @@
+cmake = import('cmake')
+
+kinsol_cmake_options = cmake.subproject_options()
+
+kinsol_cmake_options.add_cmake_defines({
+    'CMAKE_CXX_FLAGS'   : '-Wno-deprecated-declarations',
+    'CMAKE_C_FLAGS'     : '-Wno-deprecated-declarations',
+    'BUILD_SHARED_LIBS' : 'ON',
+    'BUILD_STATIC_LIBS' : 'OFF',
+})
+
+kinsol_cmake_options.add_cmake_defines({
+    'CMAKE_INSTALL_LIBDIR': get_option('libdir'),
+    'CMAKE_INSTALL_INCLUDEDIR': get_option('includedir')
+})
+
+kinsol_sp = cmake.subproject(
+    'kinsol',
+    options: kinsol_cmake_options,
+)
+
+sundials_kinsol_shared = kinsol_sp.dependency('sundials_kinsol_shared')
+
+kinsol_dep = declare_dependency(
+    dependencies: [
+        sundials_kinsol_shared,
+    ]
+)
+
diff --git a/build-config/sundials/meson.build b/build-config/sundials/meson.build
new file mode 100644
index 00000000..e5cb6298
--- /dev/null
+++ b/build-config/sundials/meson.build
@@ -0,0 +1,9 @@
+subdir('cvode')
+subdir('kinsol')
+
+sundials_dep = declare_dependency(
+    dependencies: [
+        cvode_dep,
+        kinsol_dep,
+    ],
+)
\ No newline at end of file
diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h
index 295319a6..8b2b745c 100644
--- a/src/include/gridfire/engine/views/engine_multiscale.h
+++ b/src/include/gridfire/engine/views/engine_multiscale.h
@@ -3,6 +3,10 @@
 #include "gridfire/engine/engine_abstract.h"
 #include "gridfire/engine/views/engine_view_abstract.h"
 #include "gridfire/engine/engine_graph.h"
+#include "sundials/sundials_linearsolver.h"
+#include "sundials/sundials_matrix.h"
+#include "sundials/sundials_nvector.h"
+#include "sundials/sundials_types.h"
 
 #include "unsupported/Eigen/NonLinearOptimization"
 
@@ -78,6 +82,8 @@ namespace gridfire {
          */
         explicit MultiscalePartitioningEngineView(DynamicEngine& baseEngine);
 
+        ~MultiscalePartitioningEngineView() override;
+
         /**
          * @brief Gets the list of species in the network.
          * @return A const reference to the vector of `Species` objects representing all species
@@ -611,120 +617,81 @@ namespace gridfire {
             [[nodiscard]] [[maybe_unused]] std::string toString() const;
         };
 
-        /**
-         * @brief Functor for solving QSE abundances using Eigen's nonlinear optimization.
-         *
-         * @par Purpose
-         * This struct provides the objective function (`operator()`) and its Jacobian
-         * (`df`) to Eigen's Levenberg-Marquardt solver. The goal is to find the abundances
-         * of algebraic species that make their time derivatives (`dY/dt`) equal to zero.
-         *
-         * @how
-         * - **`operator()`**: Takes a vector `v_qse` (scaled abundances of algebraic species) as input.
-         *   It constructs a full trial abundance vector `y_trial`, calls the base engine's
-         *   `calculateRHSAndEnergy`, and returns the `dY/dt` values for the algebraic species.
-         *   The solver attempts to drive this return vector to zero.
-         * - **`df`**: Computes the Jacobian of the objective function. It calls the base engine's
-         *   `generateJacobianMatrix` and extracts the sub-matrix corresponding to the algebraic
-         *   species. It applies the chain rule to account for the `asinh` scaling used on the
-         *   abundances.
-         *
-         * The abundances are scaled using `asinh` to handle the large dynamic range and ensure positivity.
-         */
-        struct EigenFunctor {
-            using InputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
-            using OutputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
-            using JacobianType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
-            enum {
-                InputsAtCompileTime = Eigen::Dynamic,
-                ValuesAtCompileTime = Eigen::Dynamic
+        class QSESolver {
+        public:
+            QSESolver(
+                const std::vector<fourdst::atomic::Species>& species,
+                const DynamicEngine& engine,
+                SUNContext sun_ctx
+            );
+
+            QSESolver(
+                const QSESolver& other
+            ) = delete;
+            QSESolver& operator=(
+                const QSESolver& other
+            ) = delete;
+
+            ~QSESolver();
+
+            fourdst::composition::Composition solve(
+                const fourdst::composition::Composition& comp,
+                double T9,
+                double rho
+            ) const;
+
+            size_t solves() const;
+
+            void log_diagnostics() const;
+        private:
+
+            static int sys_func(
+                N_Vector y,
+                N_Vector f,
+                void *user_data
+            );
+            static int sys_jac(
+                N_Vector y,
+                N_Vector fy,
+                SUNMatrix J,
+                void *user_data,
+                N_Vector tmp1,
+                N_Vector tmp2
+            );
+
+            struct UserData {
+                const DynamicEngine& engine;
+                double T9;
+                double rho;
+                fourdst::composition::Composition& comp;
+                const std::unordered_map<fourdst::atomic::Species, size_t>& qse_solve_species_index_map;
+                const std::vector<fourdst::atomic::Species>& qse_solve_species;
             };
 
-            /**
-             * @brief Pointer to the MultiscalePartitioningEngineView instance.
-             */
-            const MultiscalePartitioningEngineView& m_view;
-            /**
-             * @brief The set of species to solve for in the QSE group.
-             */
-            const std::set<fourdst::atomic::Species>& m_qse_solve_species;
-            /**
-             * @brief Initial abundances of all species in the full network.
-             */
-            const fourdst::composition::CompositionAbstract& m_initial_comp;
-            /**
-             * @brief Temperature in units of 10^9 K.
-             */
-            const double m_T9;
-            /**
-             * @brief Density in g/cm^3.
-             */
-            const double m_rho;
-            /**
-             * @brief Scaling factors for the species abundances, used to improve solver stability.
-             */
-            const Eigen::VectorXd& m_Y_scale;
+        private:
+            mutable size_t m_solves = 0;
+            size_t m_N;
+            const DynamicEngine& m_engine;
+            std::vector<fourdst::atomic::Species> m_species;
+            std::unordered_map<fourdst::atomic::Species, size_t> m_speciesMap;
 
-            /**
-             * @brief Mapping from species to their indices in the QSE solve vector.
-             */
-            const std::unordered_map<fourdst::atomic::Species, size_t> m_qse_solve_species_index_map;
+            // --- SUNDIALS / KINSOL Persistent resources
+            SUNContext m_sun_ctx   = nullptr;
+            void* m_kinsol_mem     = nullptr;
 
-            mutable std::optional<JacobianType> m_cached_jacobian = std::nullopt;
+            N_Vector m_Y           = nullptr;
+            N_Vector m_scale       = nullptr;
+            N_Vector m_f_scale     = nullptr;
+            N_Vector m_constraints = nullptr;
+            N_Vector m_func_tmpl   = nullptr;
 
-            /**
-             * @brief Constructs an EigenFunctor.
-             *
-             * @param view The MultiscalePartitioningEngineView instance.
-             * @param qse_solve_species Species to solve for in the QSE group.
-             * @param initial_comp Initial abundances of all species in the full network.
-             * @param T9 Temperature in units of 10^9 K.
-             * @param rho Density in g/cm^3.
-             * @param Y_scale Scaling factors for the species abundances.
-             * @param qse_solve_species_index_map Mapping from species to their indices in the QSE solve vector.
-             */
-            EigenFunctor(
-                const MultiscalePartitioningEngineView& view,
-                const std::set<fourdst::atomic::Species>& qse_solve_species,
-                const fourdst::composition::CompositionAbstract& initial_comp,
-                const double T9,
-                const double rho,
-                const Eigen::VectorXd& Y_scale,
-                const std::unordered_map<fourdst::atomic::Species, size_t>& qse_solve_species_index_map
-            ) :
-            m_view(view),
-            m_qse_solve_species(qse_solve_species),
-            m_initial_comp(initial_comp),
-            m_T9(T9),
-            m_rho(rho),
-            m_Y_scale(Y_scale),
-            m_qse_solve_species_index_map(qse_solve_species_index_map) {}
+            SUNMatrix m_J          = nullptr;
+            SUNLinearSolver m_LS   = nullptr;
 
-            /**
-             * @brief Gets the number of output values from the functor (size of the residual vector).
-             * @return The number of algebraic species being solved.
-             */
-            [[nodiscard]] size_t values() const { return m_qse_solve_species.size(); }
-            /**
-             * @brief Gets the number of input values to the functor (size of the variable vector).
-             * @return The number of algebraic species being solved.
-             */
-            [[nodiscard]] size_t inputs() const { return m_qse_solve_species.size(); }
+            static quill::Logger* getLogger() {
+                return LogManager::getInstance().getLogger("log");
+            }
 
-            /**
-             * @brief Evaluates the functor's residual vector `f_qse = dY_alg/dt`.
-             * @param v_qse The input vector of scaled algebraic abundances.
-             * @param f_qse The output residual vector.
-             * @return 0 on success.
-             */
-            int operator()(const InputType& v_qse, OutputType& f_qse) const;
-            /**
-             * @brief Evaluates the Jacobian of the functor, `J_qse = d(f_qse)/d(v_qse)`.
-             * @param v_qse The input vector of scaled algebraic abundances.
-             * @param J_qse The output Jacobian matrix.
-             * @return 0 on success.
-             */
-            int df(const InputType& v_qse, JacobianType& J_qse) const;
         };
 
         struct FluxValidationResult {
@@ -746,6 +713,11 @@ namespace gridfire {
          * @brief The list of identified equilibrium groups.
          */
         std::vector<QSEGroup> m_qse_groups;
+
+        /**
+         * @brief A set of solvers, one for each QSE group
+         */
+        std::vector<std::unique_ptr<QSESolver>> m_qse_solvers;
         /**
          * @brief The simplified set of species presented to the solver (the "slow" species).
          */
@@ -771,6 +743,8 @@ namespace gridfire {
 
         mutable std::unordered_map<uint64_t, fourdst::composition::Composition> m_composition_cache;
 
+        SUNContext m_sun_ctx = nullptr; // TODO: initialize and safely destroy this
+
     private:
         /**
          * @brief Partitions the network by timescale.
diff --git a/src/include/gridfire/utils/sundials.h b/src/include/gridfire/utils/sundials.h
new file mode 100644
index 00000000..a8f44cf9
--- /dev/null
+++ b/src/include/gridfire/utils/sundials.h
@@ -0,0 +1,66 @@
+#pragma once
+#include <unordered_map>
+#include <nvector/nvector_serial.h>
+
+#include "gridfire/exceptions/error_solver.h"
+#include "sundials/sundials_nvector.h"
+
+namespace gridfire::utils {
+    inline std::unordered_map<int, std::string> cvode_ret_code_map {
+        {0, "CV_SUCCESS: The solver succeeded."},
+        {1, "CV_TSTOP_RETURN: The solver reached the specified stopping time."},
+        {2, "CV_ROOT_RETURN: A root was found."},
+        {-99, "CV_WARNING: CVODE succeeded but in an unusual manner"},
+        {-1, "CV_TOO_MUCH_WORK: The solver took too many internal steps."},
+        {-2, "CV_TOO_MUCH_ACC: The solver could not satisfy the accuracy requested."},
+        {-3, "CV_ERR_FAILURE: The solver encountered a non-recoverable error."},
+        {-4, "CV_CONV_FAILURE: The solver failed to converge."},
+        {-5, "CV_LINIT_FAIL: The linear solver's initialization function failed."},
+        {-6, "CV_LSETUP_FAIL: The linear solver's setup function failed."},
+        {-7, "CV_LSOLVE_FAIL: The linear solver's solve function failed."},
+        {-8, "CV_RHSFUNC_FAIL: The right-hand side function failed in an unrecoverable manner."},
+        {-9, "CV_FIRST_RHSFUNC_ERR: The right-hand side function failed at the first call."},
+        {-10, "CV_REPTD_RHSFUNC_ERR: The right-hand side function repeatedly failed recoverable."},
+        {-11, "CV_UNREC_RHSFUNC_ERR: The right-hand side function failed unrecoverably."},
+        {-12, "CV_RTFUNC_FAIL: The rootfinding function failed in an unrecoverable manner."},
+        {-13, "CV_NLS_INIT_FAIL: The nonlinear solver's initialization function failed."},
+        {-14, "CV_NLS_SETUP_FAIL: The nonlinear solver's setup function failed."},
+        {-15, "CV_CONSTR_FAIL : The inequality constraint was violated and the solver was unable to recover."},
+        {-16, "CV_NLS_FAIL: The nonlinear solver's solve function failed."},
+        {-20, "CV_MEM_FAIL: Memory allocation failed."},
+        {-21, "CV_MEM_NULL: The CVODE memory structure is NULL."},
+        {-22, "CV_ILL_INPUT: An illegal input was detected."},
+        {-23, "CV_NO_MALLOC: The CVODE memory structure has not been allocated."},
+        {-24, "CV_BAD_K: The value of k is invalid."},
+        {-25, "CV_BAD_T: The value of t is invalid."},
+        {-26, "CV_BAD_DKY: The value of dky is invalid."},
+        {-27, "CV_TOO_CLOSE: The time points are too close together."},
+        {-28, "CV_VECTOROP_ERR: A vector operation failed."},
+        {-29, "CV_PROJ_MEM_NULL: The projection memory structure is NULL."},
+        {-30, "CV_PROJFUNC_FAIL: The projection function failed in an unrecoverable manner."},
+        {-31, "CV_REPTD_PROJFUNC_ERR: The projection function has repeated recoverable errors."}
+    };
+
+    inline void check_cvode_flag(const int flag, const std::string& func_name) {
+        if (flag < 0) {
+            if (!cvode_ret_code_map.contains(flag)) {
+                throw exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": Unknown error code: " + std::to_string(flag));
+            }
+            throw exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": " + cvode_ret_code_map.at(flag));
+        }
+    }
+
+    inline N_Vector init_sun_vector(uint64_t size, SUNContext sun_ctx) {
+        #ifdef SUNDIALS_HAVE_OPENMP
+            const N_Vector vec = N_VNew_OpenMP(size, 0, sun_ctx);
+        #elif SUNDIALS_HAVE_PTHREADS
+            const N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
+        #else
+            const N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
+        #endif
+
+        check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
+        return vec;
+    }
+
+}
diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp
index fe55141a..576e2315 100644
--- a/src/lib/engine/views/engine_multiscale.cpp
+++ b/src/lib/engine/views/engine_multiscale.cpp
@@ -1,6 +1,7 @@
 #include "gridfire/engine/views/engine_multiscale.h"
 #include "gridfire/exceptions/error_engine.h"
 #include "gridfire/engine/procedures/priming.h"
+#include "gridfire/utils/sundials.h"
 
 #include <stdexcept>
 #include <vector>
@@ -17,12 +18,15 @@
 #include "quill/LogMacros.h"
 #include "quill/Logger.h"
 
+#include "kinsol/kinsol.h"
+#include "sundials/sundials_context.h"
+#include "sunmatrix/sunmatrix_dense.h"
+#include "sunlinsol/sunlinsol_dense.h"
+
 #include  "xxhash64.h"
 #include "fourdst/composition/utils/composition_hash.h"
 
 namespace {
-    constexpr double MIN_ABS_NORM_ALWAYS_CONVERGED = 1e-30;
-
     using namespace fourdst::atomic;
 
 
@@ -159,6 +163,18 @@ namespace {
         {Eigen::LevenbergMarquardtSpace::Status::GtolTooSmall, "GtolTooSmall"},
         {Eigen::LevenbergMarquardtSpace::Status::UserAsked, "UserAsked"}
     };
+
+    void QuietErrorRouter(int line, const char *func, const char *file, const char *msg,
+                      SUNErrCode err_code, void *err_user_data, SUNContext sunctx) {
+
+        // LIST OF ERRORS TO IGNORE
+        if (err_code == KIN_LINESEARCH_NONCONV) {
+            return;
+        }
+
+        // For everything else, use the default SUNDIALS logger (or your own)
+        SUNLogErrHandlerFn(line, func, file, msg, err_code, err_user_data, sunctx);
+    }
 }
 
 namespace gridfire {
@@ -166,7 +182,23 @@ namespace gridfire {
 
     MultiscalePartitioningEngineView::MultiscalePartitioningEngineView(
         DynamicEngine& baseEngine
-    ) : m_baseEngine(baseEngine) {}
+    ) : m_baseEngine(baseEngine) {
+        const int flag = SUNContext_Create(SUN_COMM_NULL, &m_sun_ctx);
+        if (flag != 0) {
+            LOG_CRITICAL(m_logger, "Error while creating SUNContext in MultiscalePartitioningEngineView");
+            throw std::runtime_error("Error creating SUNContext in MultiscalePartitioningEngineView");
+        }
+        SUNContext_PushErrHandler(m_sun_ctx, QuietErrorRouter, nullptr);
+    }
+
+    MultiscalePartitioningEngineView::~MultiscalePartitioningEngineView() {
+        LOG_TRACE_L1(m_logger, "Cleaning up MultiscalePartitioningEngineView...");
+        m_qse_solvers.clear();
+        if (m_sun_ctx) {
+            SUNContext_Free(&m_sun_ctx);
+            m_sun_ctx = nullptr;
+        }
+    }
 
     const std::vector<Species> & MultiscalePartitioningEngineView::getNetworkSpecies() const {
         return m_baseEngine.getNetworkSpecies();
@@ -598,6 +630,7 @@ namespace gridfire {
         LOG_TRACE_L1(m_logger, "Partitioning network...");
         LOG_TRACE_L1(m_logger, "Clearing previous state...");
         m_qse_groups.clear();
+        m_qse_solvers.clear();
         m_dynamic_species.clear();
         m_algebraic_species.clear();
         m_composition_cache.clear(); // We need to clear the cache now cause the same comp, temp, and density may result in a different value
@@ -790,7 +823,18 @@ namespace gridfire {
                 m_dynamic_species.push_back(species);
             }
         }
-        return getNormalizedEquilibratedComposition(comp, T9, rho);
+
+        for (const auto& group : m_qse_groups) {
+            std::vector<Species> groupAlgebraicSpecies;
+            for (const auto& species : group.algebraic_species) {
+                groupAlgebraicSpecies.push_back(species);
+            }
+            m_qse_solvers.push_back(std::make_unique<QSESolver>(groupAlgebraicSpecies, m_baseEngine, m_sun_ctx));
+        }
+
+        fourdst::composition::Composition result = getNormalizedEquilibratedComposition(comp, T9, rho);
+
+        return result;
     }
 
     void MultiscalePartitioningEngineView::exportToDot(
@@ -1502,98 +1546,18 @@ namespace gridfire {
         const double T9,
         const double rho
     ) const {
-        LOG_TRACE_L1(m_logger, "Solving for QSE abundances...");
-        LOG_TRACE_L2(m_logger, "Composition before QSE solving: {}", [&comp]() -> std::string {
-            std::stringstream ss;
-            size_t i = 0;
-            for (const auto& [sp, y] : comp) {
-                ss << std::format("{}: {}", sp.name(), y);
-                if (i < comp.size() - 1) {
-                    ss << ", ";
-                }
-                i++;
-            }
-            return ss.str();
-        }());
+        LOG_TRACE_L2(m_logger, "Solving for QSE abundances...");
 
         fourdst::composition::Composition outputComposition(comp);
 
-        for (const auto&[is_in_equilibrium, algebraic_species, seed_species, mean_timescale] : m_qse_groups) {
-            LOG_TRACE_L2(m_logger, "Working on QSE group with algebraic species: {}",
-                [&]() -> std::string {
-                    std::stringstream ss;
-                    size_t count = 0;
-                    for (const auto& species: algebraic_species) {
-                        ss << species.name();
-                        if (count < algebraic_species.size() - 1) {
-                            ss << ", ";
-                        }
-                        count++;
-                    }
-                    return ss.str();
-                }());
-            if (!is_in_equilibrium || (algebraic_species.empty() && seed_species.empty())) {
-                continue;
-            }
-
-            Eigen::VectorXd Y_scale(algebraic_species.size());
-            Eigen::VectorXd v_initial(algebraic_species.size());
-            long i = 0;
-            std::unordered_map<Species, size_t> species_to_index_map;
-            for (const auto& species : algebraic_species) {
-                constexpr double abundance_floor = 1.0e-100;
-                const double initial_abundance = comp.getMolarAbundance(species);
-                const double Y = std::max(initial_abundance, abundance_floor);
-                v_initial(i) = std::log(Y);
-                species_to_index_map.emplace(species, i);
-                LOG_TRACE_L2(m_logger, "For species {} initial molar abundance is {}, log scaled to {}. Species placed at index {}.", species.name(), Y, v_initial(i), i);
-                i++;
-            }
-
-            LOG_TRACE_L2(m_logger, "Setting up Eigen Levenberg-Marquardt solver for QSE group...");
-            EigenFunctor functor(*this, algebraic_species, comp, T9, rho, Y_scale, species_to_index_map);
-            Eigen::LevenbergMarquardt lm(functor);
-            lm.parameters.ftol = 1.0e-10;
-            lm.parameters.xtol = 1.0e-10;
-
-            LOG_TRACE_L2(m_logger, "Minimizing functor...");
-            Eigen::LevenbergMarquardtSpace::Status status = lm.minimize(v_initial);
-            LOG_TRACE_L2(m_logger, "Minimizing functor status: {}", lm_status_map.at(status));
-
-            if (status <= 0 || status > 4) {
-                std::stringstream msg;
-                msg << "While working on QSE group with algebraic species: ";
-                size_t count = 0;
-                for (const auto& species: algebraic_species) {
-                    msg << species;
-                    if (count < algebraic_species.size() - 1) {
-                        msg << ", ";
-                    }
-                    count++;
-                }
-                msg << " the QSE solver failed to converge with status: " << lm_status_map.at(status) << " (No. " << status << ").";
-                LOG_ERROR(m_logger, "{}", msg.str());
-                throw std::runtime_error(msg.str());
-            }
-            LOG_TRACE_L1(m_logger, "QSE Group minimization succeeded with status: {}", lm_status_map.at(status));
-            Eigen::VectorXd Y_final_qse = v_initial.array().exp(); // Convert back to physical abundances using exponential scaling
-            i = 0;
-            for (const auto& species: algebraic_species) {
-                LOG_TRACE_L1(
-                    m_logger,
-                    "During QSE solving species {} started with a molar abundance of {} and ended with an abundance of {}.",
-                    species.name(),
-                    comp.getMolarAbundance(species),
-                    Y_final_qse(i)
-                );
-                // double Xi = Y_final_qse(i) * species.mass(); // Convert from molar abundance to mass fraction
-                if (!outputComposition.contains (species)) {
-                    outputComposition.registerSpecies(species);
-                }
-                outputComposition.setMolarAbundance(species, Y_final_qse(i));
-                i++;
+        for (const auto& [group, solver]: std::views::zip(m_qse_groups, m_qse_solvers)) {
+            const fourdst::composition::Composition groupResult = solver->solve(outputComposition, T9, rho);
+            for (const auto& [sp, y] : groupResult) {
+                outputComposition.setMolarAbundance(sp, y);
             }
+            solver->log_diagnostics();
         }
+        LOG_TRACE_L2(m_logger, "Done solving for QSE abundances!");
         return outputComposition;
     }
 
@@ -1795,121 +1759,249 @@ namespace gridfire {
         return candidate_groups;
     }
 
-
-     //////////////////////////////////////
-    /// Eigen Functor Member Functions ///
-   /////////////////////////////////////
+     //////////////////////////////////
+    /// QSESolver Member Functions ///
+    //////////////////////////////////
 
 
-    int MultiscalePartitioningEngineView::EigenFunctor::operator()(const InputType &v_qse, OutputType &f_qse) const {
-        fourdst::composition::Composition comp_trial(m_initial_comp.getRegisteredSpecies());
-        for (const auto& [sp, y] : m_initial_comp) {
-            comp_trial.setMolarAbundance(sp, y);
-        }
-        Eigen::VectorXd y_qse = v_qse.array().exp(); // Convert to physical abundances using exponential scaling
+    MultiscalePartitioningEngineView::QSESolver::QSESolver(
+        const std::vector<fourdst::atomic::Species>& species,
+        const DynamicEngine& engine,
+        const SUNContext sun_ctx
+    ) :
+    m_N(species.size()),
+    m_engine(engine),
+    m_species(species),
+    m_sun_ctx(sun_ctx) {
+        m_Y = utils::init_sun_vector(m_N, m_sun_ctx);
+        m_scale = N_VClone(m_Y);
+        m_f_scale = N_VClone(m_Y);
+        m_constraints = N_VClone(m_Y);
+        m_func_tmpl = N_VClone(m_Y);
 
-        for (const auto& species: m_qse_solve_species) {
-            auto index = static_cast<long>(m_qse_solve_species_index_map.at(species));
-            comp_trial.setMolarAbundance(species, y_qse[index]);
+        if (!m_Y || !m_scale || !m_constraints || !m_func_tmpl) {
+            LOG_CRITICAL(getLogger(), "Failed to allocate SUNVectors for QSE solver.");
+            throw std::runtime_error("Failed to allocate SUNVectors for QSE solver.");
         }
 
-        const auto result = m_view.getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho);
-        if (!result) {
-            throw exceptions::StaleEngineError("Failed to calculate RHS and energy due to stale engine state");
+        for (size_t i = 0; i < m_N; ++i) {
+            m_speciesMap[m_species[i]] = i;
         }
-        const auto&[dydt, nuclearEnergyGenerationRate, _] = result.value();
-        f_qse.resize(static_cast<long>(m_qse_solve_species.size()));
-        long i = 0;
-        // TODO: make sure that just counting up i is a valid approach, this is a possible place an indexing bug may have crept in
-        for (const auto& species: m_qse_solve_species) {
-            const double dydt_i = dydt.at(species);
-            f_qse(i) = dydt_i/y_qse(i); // We square the residuals to improve numerical stability in the solver
-            i++;
-        }
-        LOG_TRACE_L2(
-            m_view.m_logger,
-            "Functor evaluation at T9 = {}, rho = {}, y_qse (v_qse) = <{}> complete. ||f|| = {}",
-            m_T9,
-            m_rho,
-            [&]() -> std::string {
-                std::stringstream ss;
-                for (long j = 0; j < y_qse.size(); ++j) {
-                    ss << y_qse(j);
-                    ss << "(" << v_qse(j) << ")";
-                    if (j < y_qse.size() - 1) {
-                        ss << ", ";
-                    }
-                }
-                return ss.str();
-            }(),
-            f_qse.norm()
-        );
-        LOG_TRACE_L3(
-            m_view.m_logger,
-            "{}",
-            [&]() -> std::string {
-                std::stringstream ss;
-                const std::vector species(m_qse_solve_species.begin(), m_qse_solve_species.end());
-                for (long j = 0; j < f_qse.size(); ++j) {
-                    ss << "Residual for species " << species.at(j).name() << " f(" << j << ") = " << f_qse(j) << "\n";
-                }
-                return ss.str();
-            }()
-        );
-        return 0;
+
+        N_VConst(1.0, m_constraints);
+        m_kinsol_mem = KINCreate(m_sun_ctx);
+
+        utils::check_cvode_flag(m_kinsol_mem ? 0 : -1, "KINCreate");
+        utils::check_cvode_flag(KINInit(m_kinsol_mem, sys_func, m_func_tmpl), "KINInit");
+        utils::check_cvode_flag(KINSetConstraints(m_kinsol_mem, m_constraints), "KINSetConstraints");
+
+        m_J = SUNDenseMatrix(static_cast<sunindextype>(m_N), static_cast<sunindextype>(m_N), m_sun_ctx);
+        utils::check_cvode_flag(m_J ? 0 : -1, "SUNDenseMatrix");
+
+        m_LS = SUNLinSol_Dense(m_Y, m_J, m_sun_ctx);
+        utils::check_cvode_flag(m_LS ? 0 : -1, "SUNLinSol_Dense");
+
+        utils::check_cvode_flag(KINSetLinearSolver(m_kinsol_mem, m_LS, m_J), "KINSetLinearSolver");
+        utils::check_cvode_flag(KINSetJacFn(m_kinsol_mem, sys_jac), "KINSetJacFn");
+
+        utils::check_cvode_flag(KINSetMaxSetupCalls(m_kinsol_mem, 20), "KINSetMaxSetupCalls");
+
+        utils::check_cvode_flag(KINSetFuncNormTol(m_kinsol_mem, 1e-6), "KINSetFuncNormTol");
+        utils::check_cvode_flag(KINSetNumMaxIters(m_kinsol_mem, 200), "KINSetNumMaxIters");
+
+        // We want to effectively disable this since enormous changes in order of magnitude are realistic for this problem.
+        utils::check_cvode_flag(KINSetMaxNewtonStep(m_kinsol_mem, 200), "KINSetMaxNewtonStep");
+
     }
 
-    int MultiscalePartitioningEngineView::EigenFunctor::df(const InputType &v_qse, JacobianType &J_qse) const {
-        fourdst::composition::Composition comp_trial(m_initial_comp.getRegisteredSpecies());
-        for (const auto& [sp, y] : m_initial_comp) {
-            comp_trial.setMolarAbundance(sp, y);
+    MultiscalePartitioningEngineView::QSESolver::~QSESolver() {
+        if (m_Y) {
+            N_VDestroy(m_Y);
+            m_Y = nullptr;
         }
-        Eigen::VectorXd y_qse = v_qse.array().exp(); // Convert to physical abundances using exponential scaling
+        if (m_scale) {
+            N_VDestroy(m_scale);
+            m_scale = nullptr;
+        }
+        if (m_f_scale) {
+            N_VDestroy(m_f_scale);
+            m_f_scale = nullptr;
+        }
+        if (m_constraints) {
+            N_VDestroy(m_constraints);
+            m_constraints = nullptr;
+        }
+        if (m_func_tmpl) {
+            N_VDestroy(m_func_tmpl);
+            m_func_tmpl = nullptr;
+        }
+        if (m_kinsol_mem) {
+            KINFree(&m_kinsol_mem);
+            m_kinsol_mem = nullptr;
+        }
+        if (m_J) {
+            SUNMatDestroy(m_J);
+            m_J = nullptr;
+        }
+        if (m_LS) {
+            SUNLinSolFree(m_LS);
+            m_LS = nullptr;
+        }
+    }
 
-        for (const auto& species: m_qse_solve_species) {
-            const double molarAbundance = y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))];
-            comp_trial.setMolarAbundance(species, molarAbundance);
+    fourdst::composition::Composition MultiscalePartitioningEngineView::QSESolver::solve(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho
+    ) const {
+
+        fourdst::composition::Composition result = comp;
+        UserData data {
+            m_engine,
+            T9,
+            rho,
+            result,
+            m_speciesMap,
+            m_species
+        };
+
+        utils::check_cvode_flag(KINSetUserData(m_kinsol_mem, &data), "KINSetUserData");
+
+        sunrealtype* y_data = N_VGetArrayPointer(m_Y);
+        sunrealtype* scale_data = N_VGetArrayPointer(m_scale);
+
+        // It is more cache optimized to do a standard as opposed to range based for-loop here
+        for (size_t i = 0; i < m_N; ++i) {
+            const auto& species = m_species[i];
+            double Y = result.getMolarAbundance(species);
+
+            constexpr double abundance_floor = 1.0e-100;
+            Y = std::max(abundance_floor, Y);
+            y_data[i] = Y;
+            scale_data[i] = 1.0;
         }
 
-        std::vector<Species> qse_species_vector(m_qse_solve_species.begin(), m_qse_solve_species.end());
-        NetworkJacobian jac = m_view.getBaseEngine().generateJacobianMatrix(comp_trial, m_T9, m_rho, qse_species_vector);
-        const auto result = m_view.getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho);
+        auto initial_rhs = m_engine.calculateRHSAndEnergy(result, T9, rho);
+        if (!initial_rhs) {
+            throw std::runtime_error("In QSE solver failed to calculate initial RHS");
+        }
+
+        sunrealtype* f_scale_data = N_VGetArrayPointer(m_f_scale);
+        for (size_t i = 0; i < m_N; ++i) {
+            const auto& species = m_species[i];
+            double dydt = std::abs(initial_rhs.value().dydt.at(species));
+            f_scale_data[i] = 1.0 / (dydt + 1e-15);
+        }
+
+        if (m_solves > 0) {
+            // After the initial solution we want to allow kinsol to reuse its state
+            utils::check_cvode_flag(KINSetNoInitSetup(m_kinsol_mem, SUNTRUE), "KINSetNoInitSetup");
+        } else {
+            utils::check_cvode_flag(KINSetNoInitSetup(m_kinsol_mem, SUNFALSE), "KINSetNoInitSetup");
+        }
+
+        const int flag = KINSol(m_kinsol_mem, m_Y, KIN_LINESEARCH, m_scale, m_f_scale);
+        if (flag < 0) {
+            LOG_WARNING(getLogger(), "KINSol failed to converge while solving QSE abundances with flag {}.", utils::cvode_ret_code_map.at(flag));
+            return comp;
+        }
+
+        for (size_t i = 0; i < m_N; ++i) {
+            const auto& species = m_species[i];
+            result.setMolarAbundance(species, y_data[i]);
+        }
+
+        m_solves++;
+        return result;
+    }
+
+    size_t MultiscalePartitioningEngineView::QSESolver::solves() const {
+        return m_solves;
+    }
+
+    void MultiscalePartitioningEngineView::QSESolver::log_diagnostics() const {
+        long int nni, nfe, nje;
+
+        int flag = KINGetNumNonlinSolvIters(m_kinsol_mem, &nni);
+
+        flag = KINGetNumFuncEvals(m_kinsol_mem, &nfe);
+
+        flag = KINGetNumJacEvals(m_kinsol_mem, &nje);
+
+        LOG_INFO(getLogger(),
+            "QSE Stats | Iters: {} | RHS Evals: {} | Jac Evals: {} | Ratio (J/I): {:.2f}",
+            nni, nfe, nje, static_cast<double>(nje) / static_cast<double>(nni)
+        );
+        getLogger()->flush_log(true);
+    }
+
+
+    int MultiscalePartitioningEngineView::QSESolver::sys_func(
+        const N_Vector y,
+        const N_Vector f,
+        void *user_data
+    ) {
+        const auto* data = static_cast<UserData*>(user_data);
+
+        const sunrealtype* y_data = N_VGetArrayPointer(y);
+        sunrealtype* f_data = N_VGetArrayPointer(f);
+
+        const auto& map = data->qse_solve_species_index_map;
+        for (size_t index = 0; index < data->qse_solve_species.size(); ++index) {
+            const auto& species = data->qse_solve_species[index];
+            data->comp.setMolarAbundance(species, y_data[index]);
+        }
+
+        const auto result = data->engine.calculateRHSAndEnergy(data->comp, data->T9, data->rho);
+
         if (!result) {
-            throw exceptions::StaleEngineError("Failed to calculate RHS and energy due to stale engine state");
-        }
-        const auto&[dydt, nuclearEnergyGenerationRate, _] = result.value();
-
-        const long N = static_cast<long>(m_qse_solve_species.size());
-        J_qse.resize(N, N);
-        long rowID = 0;
-        for (const auto& rowSpecies : m_qse_solve_species) {
-            long colID = 0;
-            for (const auto& colSpecies: m_qse_solve_species) {
-                J_qse(rowID, colID) = jac(rowSpecies, colSpecies);
-                colID += 1;
-                LOG_TRACE_L3(m_view.m_logger, "Jacobian[{}, {}] (d(dY({}))/dY({})) = {}", rowID, colID - 1, rowSpecies.name(), colSpecies.name(), J_qse(rowID, colID - 1));
-            }
-            rowID += 1;
+            return 1; // Potentially recoverable error
         }
 
-        for (long i = 0; i < J_qse.rows(); ++i) {
-            for (long j = 0; j < J_qse.cols(); ++j) {
-                double on_diag_correction = 0.0;
-                if (i == j) {
-                    auto rowSpecies = *(std::next(m_qse_solve_species.begin(), i));
-                    const double Fi = dydt.at(rowSpecies);
-                    on_diag_correction = Fi / y_qse(i);
-                }
-                J_qse(i, j) = y_qse(j) * (J_qse(i, j) - on_diag_correction) / y_qse(i); // Apply chain rule J'(i,j) = y_j * (J(i,j) - δ_ij(F_i/Y_i)) / Y_i
-            }
-        }
+        const auto& dydt = result.value().dydt;
 
-        m_cached_jacobian = J_qse; // Cache the computed Jacobian for future use
+        for (const auto& [species, index] : map) {
+            f_data[index] = dydt.at(species);
+        }
 
         return 0; // Success
     }
 
 
+    int MultiscalePartitioningEngineView::QSESolver::sys_jac(
+        const N_Vector y,
+        N_Vector fy,
+        SUNMatrix J,
+        void *user_data,
+        N_Vector tmp1,
+        N_Vector tmp2
+    ) {
+        const auto* data = static_cast<UserData*>(user_data);
+
+        const sunrealtype* y_data = N_VGetArrayPointer(y);
+        const auto& map = data->qse_solve_species_index_map;
+        for (const auto& [species, index] : map) {
+            data->comp.setMolarAbundance(species, y_data[index]);
+        }
+
+        const NetworkJacobian jac = data->engine.generateJacobianMatrix(
+            data->comp,
+            data->T9,
+            data->rho,
+            data->qse_solve_species
+        );
+
+        sunrealtype* J_data = SUNDenseMatrix_Data(J);
+        const sunindextype N = SUNDenseMatrix_Columns(J);
+
+        for (const auto& [col_species, col_idx] : map) {
+            for (const auto& [row_species, row_idx] : map) {
+                J_data[col_idx * N + row_idx] = jac(row_species, col_species);
+            }
+        }
+
+        return 0;
+    }
 
       /////////////////////////////////
      /// QSEGroup Member Functions ///
diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp
index ff0d6093..7ec7e16c 100644
--- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp
+++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp
@@ -22,65 +22,8 @@
 #include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h"
 #include "gridfire/trigger/procedures/trigger_pprint.h"
 #include "gridfire/exceptions/error_solver.h"
-#include "gridfire/utils/logging.h"
+#include "gridfire/utils/sundials.h"
 
-namespace {
-
-    std::unordered_map<int, std::string> cvode_ret_code_map {
-        {0, "CV_SUCCESS: The solver succeeded."},
-        {1, "CV_TSTOP_RETURN: The solver reached the specified stopping time."},
-        {2, "CV_ROOT_RETURN: A root was found."},
-        {-99, "CV_WARNING: CVODE succeeded but in an unusual manner"},
-        {-1, "CV_TOO_MUCH_WORK: The solver took too many internal steps."},
-        {-2, "CV_TOO_MUCH_ACC: The solver could not satisfy the accuracy requested."},
-        {-3, "CV_ERR_FAILURE: The solver encountered a non-recoverable error."},
-        {-4, "CV_CONV_FAILURE: The solver failed to converge."},
-        {-5, "CV_LINIT_FAIL: The linear solver's initialization function failed."},
-        {-6, "CV_LSETUP_FAIL: The linear solver's setup function failed."},
-        {-7, "CV_LSOLVE_FAIL: The linear solver's solve function failed."},
-        {-8, "CV_RHSFUNC_FAIL: The right-hand side function failed in an unrecoverable manner."},
-        {-9, "CV_FIRST_RHSFUNC_ERR: The right-hand side function failed at the first call."},
-        {-10, "CV_REPTD_RHSFUNC_ERR: The right-hand side function repeatedly failed recoverable."},
-        {-11, "CV_UNREC_RHSFUNC_ERR: The right-hand side function failed unrecoverably."},
-        {-12, "CV_RTFUNC_FAIL: The rootfinding function failed in an unrecoverable manner."},
-        {-13, "CV_NLS_INIT_FAIL: The nonlinear solver's initialization function failed."},
-        {-14, "CV_NLS_SETUP_FAIL: The nonlinear solver's setup function failed."},
-        {-15, "CV_CONSTR_FAIL : The inequality constraint was violated and the solver was unable to recover."},
-        {-16, "CV_NLS_FAIL: The nonlinear solver's solve function failed."},
-        {-20, "CV_MEM_FAIL: Memory allocation failed."},
-        {-21, "CV_MEM_NULL: The CVODE memory structure is NULL."},
-        {-22, "CV_ILL_INPUT: An illegal input was detected."},
-        {-23, "CV_NO_MALLOC: The CVODE memory structure has not been allocated."},
-        {-24, "CV_BAD_K: The value of k is invalid."},
-        {-25, "CV_BAD_T: The value of t is invalid."},
-        {-26, "CV_BAD_DKY: The value of dky is invalid."},
-        {-27, "CV_TOO_CLOSE: The time points are too close together."},
-        {-28, "CV_VECTOROP_ERR: A vector operation failed."},
-        {-29, "CV_PROJ_MEM_NULL: The projection memory structure is NULL."},
-        {-30, "CV_PROJFUNC_FAIL: The projection function failed in an unrecoverable manner."},
-        {-31, "CV_REPTD_PROJFUNC_ERR: The projection function has repeated recoverable errors."}
-    };
-    void check_cvode_flag(const int flag, const std::string& func_name) {
-        if (flag < 0) {
-            if (!cvode_ret_code_map.contains(flag)) {
-                throw gridfire::exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": Unknown error code: " + std::to_string(flag));
-            }
-            throw gridfire::exceptions::CVODESolverFailureError("CVODE error in " + func_name + ": " + cvode_ret_code_map.at(flag));
-        }
-    }
-
-    N_Vector init_sun_vector(uint64_t size, SUNContext sun_ctx) {
-#ifdef SUNDIALS_HAVE_OPENMP
-        const N_Vector vec = N_VNew_OpenMP(size, 0, sun_ctx);
-#elif SUNDIALS_HAVE_PTHREADS
-        const N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
-#else
-        const N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
-#endif
-        check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
-        return vec;
-    }
-}
 
 namespace gridfire::solver {
 
@@ -175,7 +118,7 @@ namespace gridfire::solver {
         LOG_TRACE_L1(m_logger, "Number of species: {} ({} independent variables)", numSpecies, N);
         LOG_TRACE_L1(m_logger, "Initializing CVODE resources");
         m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx);
-        check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
+        utils::check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
 
         initialize_cvode_integration_resources(N, numSpecies, 0.0, equilibratedComposition, absTol, relTol, 0.0);
 
@@ -206,11 +149,11 @@ namespace gridfire::solver {
             user_data.networkSpecies = &m_engine.getNetworkSpecies();
             user_data.captured_exception.reset();
 
-            check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData");
+            utils::check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData");
 
             LOG_TRACE_L2(m_logger, "Taking one CVODE step...");
             int flag = CVode(m_cvode_mem, netIn.tMax, m_Y, &current_time, CV_ONE_STEP);
-            LOG_TRACE_L2(m_logger, "CVODE step complete. Current time: {}, step status: {}", current_time, cvode_ret_code_map.at(flag));
+            LOG_TRACE_L2(m_logger, "CVODE step complete. Current time: {}, step status: {}", current_time, utils::cvode_ret_code_map.at(flag));
 
             if (user_data.captured_exception){
                 std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception));
@@ -220,7 +163,7 @@ namespace gridfire::solver {
             // TODO: Come up with some way to save these to a file rather than spamming stdout. JSON maybe? OPAT?
             // log_step_diagnostics(user_data, true, false, false);
             // exit(0);
-            check_cvode_flag(flag, "CVode");
+            utils::check_cvode_flag(flag, "CVode");
 
             long int n_steps;
             double last_step_size;
@@ -449,11 +392,11 @@ namespace gridfire::solver {
                 cleanup_cvode_resources(true);
 
                 m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx);
-                check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
+                utils::check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate");
 
                 initialize_cvode_integration_resources(N, numSpecies, current_time, currentComposition, absTol, relTol, accumulated_energy);
 
-                check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit");
+                utils::check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit");
             }
 
         }
@@ -487,7 +430,7 @@ namespace gridfire::solver {
         NetOut netOut;
         netOut.composition = outputComposition;
         netOut.energy = accumulated_energy;
-        check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast<long int *>(&netOut.num_steps)), "CVodeGetNumSteps");
+        utils::check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast<long int *>(&netOut.num_steps)), "CVodeGetNumSteps");
 
         LOG_TRACE_L2(m_logger, "generating final nuclear energy generation rate derivatives...");
         auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives(
@@ -688,7 +631,7 @@ namespace gridfire::solver {
         std::vector<double> y_vec(y_data, y_data + numSpecies);
         fourdst::composition::Composition composition(m_engine.getNetworkSpecies(), y_vec);
 
-        LOG_TRACE_L2(m_logger, "Calculating RHS at time {} with {} species in composition (mean molecular mass: {})", t, composition.size(), composition.getMeanParticleMass());
+        LOG_TRACE_L2(m_logger, "Calculating RHS at time {} with {} species in composition", t, composition.size());
         const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho);
         if (!result) {
             LOG_WARNING(m_logger, "StaleEngineTrigger thrown during RHS calculation at time {}", t);
@@ -732,7 +675,7 @@ namespace gridfire::solver {
         LOG_TRACE_L2(m_logger, "Initializing CVODE integration resources with N: {}, current_time: {}, absTol: {}, relTol: {}", N, current_time, absTol, relTol);
         cleanup_cvode_resources(false); // Cleanup any existing resources before initializing new ones
 
-        m_Y = init_sun_vector(N, m_sun_ctx);
+        m_Y = utils::init_sun_vector(N, m_sun_ctx);
         m_YErr = N_VClone(m_Y);
 
         sunrealtype *y_data = N_VGetArrayPointer(m_Y);
@@ -747,8 +690,8 @@ namespace gridfire::solver {
         y_data[numSpecies] = accumulatedEnergy; // Specific energy rate, initialized to zero
 
 
-        check_cvode_flag(CVodeInit(m_cvode_mem, cvode_rhs_wrapper, current_time, m_Y), "CVodeInit");
-        check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances");
+        utils::check_cvode_flag(CVodeInit(m_cvode_mem, cvode_rhs_wrapper, current_time, m_Y), "CVodeInit");
+        utils::check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances");
 
         // Constraints
         // We constrain the solution vector using CVODE's built in constraint flags as outlines on page 53 of the CVODE manual
@@ -768,17 +711,17 @@ namespace gridfire::solver {
         }
         N_VConst(1.0, m_constraints); // Set all constraints to >= 0 (note this is where the flag values are set)
 
-        check_cvode_flag(CVodeSetConstraints(m_cvode_mem, m_constraints), "CVodeSetConstraints");
+        utils::check_cvode_flag(CVodeSetConstraints(m_cvode_mem, m_constraints), "CVodeSetConstraints");
 
-        check_cvode_flag(CVodeSetMaxStep(m_cvode_mem, 1.0e20), "CVodeSetMaxStep");
+        utils::check_cvode_flag(CVodeSetMaxStep(m_cvode_mem, 1.0e20), "CVodeSetMaxStep");
 
         m_J = SUNDenseMatrix(static_cast<sunindextype>(N), static_cast<sunindextype>(N), m_sun_ctx);
-        check_cvode_flag(m_J == nullptr ? -1 : 0, "SUNDenseMatrix");
+        utils::check_cvode_flag(m_J == nullptr ? -1 : 0, "SUNDenseMatrix");
         m_LS = SUNLinSol_Dense(m_Y, m_J, m_sun_ctx);
-        check_cvode_flag(m_LS == nullptr ? -1 : 0, "SUNLinSol_Dense");
+        utils::check_cvode_flag(m_LS == nullptr ? -1 : 0, "SUNLinSol_Dense");
 
-        check_cvode_flag(CVodeSetLinearSolver(m_cvode_mem, m_LS, m_J), "CVodeSetLinearSolver");
-        check_cvode_flag(CVodeSetJacFn(m_cvode_mem, cvode_jac_wrapper), "CVodeSetJacFn");
+        utils::check_cvode_flag(CVodeSetLinearSolver(m_cvode_mem, m_LS, m_J), "CVodeSetLinearSolver");
+        utils::check_cvode_flag(CVodeSetJacFn(m_cvode_mem, cvode_jac_wrapper), "CVodeSetJacFn");
         LOG_TRACE_L2(m_logger, "CVODE solver initialized");
     }
 
@@ -814,10 +757,10 @@ namespace gridfire::solver {
         sunrealtype hlast, hcur, tcur;
         int qlast;
 
-        check_cvode_flag(CVodeGetLastStep(m_cvode_mem, &hlast), "CVodeGetLastStep");
-        check_cvode_flag(CVodeGetCurrentStep(m_cvode_mem, &hcur), "CVodeGetCurrentStep");
-        check_cvode_flag(CVodeGetLastOrder(m_cvode_mem, &qlast), "CVodeGetLastOrder");
-        check_cvode_flag(CVodeGetCurrentTime(m_cvode_mem, &tcur), "CVodeGetCurrentTime");
+        utils::check_cvode_flag(CVodeGetLastStep(m_cvode_mem, &hlast), "CVodeGetLastStep");
+        utils::check_cvode_flag(CVodeGetCurrentStep(m_cvode_mem, &hcur), "CVodeGetCurrentStep");
+        utils::check_cvode_flag(CVodeGetLastOrder(m_cvode_mem, &qlast), "CVodeGetLastOrder");
+        utils::check_cvode_flag(CVodeGetCurrentTime(m_cvode_mem, &tcur), "CVodeGetCurrentTime");
 
         {
             std::vector<std::string> labels = {"Current Time (tcur)", "Last Step (hlast)", "Current Step (hcur)", "Last Order (qlast)"};
@@ -834,13 +777,13 @@ namespace gridfire::solver {
         // These are the CRITICAL counters for diagnosing your problem
         long int nsteps, nfevals, nlinsetups, netfails, nniters, nconvfails, nsetfails;
 
-        check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, &nsteps), "CVodeGetNumSteps");
-        check_cvode_flag(CVodeGetNumRhsEvals(m_cvode_mem, &nfevals), "CVodeGetNumRhsEvals");
-        check_cvode_flag(CVodeGetNumLinSolvSetups(m_cvode_mem, &nlinsetups), "CVodeGetNumLinSolvSetups");
-        check_cvode_flag(CVodeGetNumErrTestFails(m_cvode_mem, &netfails), "CVodeGetNumErrTestFails");
-        check_cvode_flag(CVodeGetNumNonlinSolvIters(m_cvode_mem, &nniters), "CVodeGetNumNonlinSolvIters");
-        check_cvode_flag(CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nconvfails), "CVodeGetNumNonlinSolvConvFails");
-        check_cvode_flag(CVodeGetNumLinConvFails(m_cvode_mem, &nsetfails), "CVodeGetNumLinConvFails");
+        utils::check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, &nsteps), "CVodeGetNumSteps");
+        utils::check_cvode_flag(CVodeGetNumRhsEvals(m_cvode_mem, &nfevals), "CVodeGetNumRhsEvals");
+        utils::check_cvode_flag(CVodeGetNumLinSolvSetups(m_cvode_mem, &nlinsetups), "CVodeGetNumLinSolvSetups");
+        utils::check_cvode_flag(CVodeGetNumErrTestFails(m_cvode_mem, &netfails), "CVodeGetNumErrTestFails");
+        utils::check_cvode_flag(CVodeGetNumNonlinSolvIters(m_cvode_mem, &nniters), "CVodeGetNumNonlinSolvIters");
+        utils::check_cvode_flag(CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nconvfails), "CVodeGetNumNonlinSolvConvFails");
+        utils::check_cvode_flag(CVodeGetNumLinConvFails(m_cvode_mem, &nsetfails), "CVodeGetNumLinConvFails");
 
 
         {
@@ -864,7 +807,7 @@ namespace gridfire::solver {
         }
 
         // --- 3. Get Estimated Local Errors (Your Original Logic) ---
-        check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors");
+        utils::check_cvode_flag(CVodeGetEstLocalErrors(m_cvode_mem, m_YErr), "CVodeGetEstLocalErrors");
 
         sunrealtype *y_data = N_VGetArrayPointer(m_Y);
         sunrealtype *y_err_data = N_VGetArrayPointer(m_YErr);
diff --git a/subprojects/cvode.wrap b/subprojects/cvode.wrap
index 99ef269e..bbbee485 100644
--- a/subprojects/cvode.wrap
+++ b/subprojects/cvode.wrap
@@ -1,6 +1,6 @@
 [wrap-file]
-directory = cvode-7.3.0
+directory = cvode-7.5.0
 
-source_url = https://github.com/LLNL/sundials/releases/download/v7.3.0/cvode-7.3.0.tar.gz
-source_filename = cvode-7.3.0.tar.gz
-source_hash = 8b15a646882f2414b1915cad4d53136717a077539e7cfc480f2002c5898ae568
+source_url = https://github.com/LLNL/sundials/releases/download/v7.5.0/cvode-7.5.0.tar.gz
+source_filename = cvode-7.5.0.tar.gz
+source_hash = e588d1c569c2fe92403c6b1a8bcce973624bf6ead4a6a302ab7cc3c4f39523bd
diff --git a/subprojects/kinsol.wrap b/subprojects/kinsol.wrap
new file mode 100644
index 00000000..abfff325
--- /dev/null
+++ b/subprojects/kinsol.wrap
@@ -0,0 +1,6 @@
+[wrap-file]
+directory = kinsol-7.5.0
+
+source_url = https://github.com/LLNL/sundials/releases/download/v7.5.0/kinsol-7.5.0.tar.gz
+source_filename = kinsol-7.5.0.tar.gz
+source_hash = ce712ea29aaec537a6cd234d1a6bea024592ef54a26a5d457dda2320f0c49d07
\ No newline at end of file
diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp
index 1b340381..85490c0f 100644
--- a/tests/graphnet_sandbox/main.cpp
+++ b/tests/graphnet_sandbox/main.cpp
@@ -66,7 +66,7 @@ gridfire::NetIn init(const double temp) {
     netIn.density = 1.5e2;
     netIn.energy = 0;
 
-    netIn.tMax = 3e17;
+    netIn.tMax = 1e17;
     netIn.dt0 = 1e-12;
 
     return netIn;