diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp
index f25ebffe9..5c3fa7dd6 100644
--- a/include/micm/solver/rosenbrock.hpp
+++ b/include/micm/solver/rosenbrock.hpp
@@ -204,7 +204,6 @@ namespace micm
       double final_time_{};
     };
 
-    System system_;
     std::vector<Process> processes_;
     RosenbrockSolverParameters parameters_;
     StateParameters state_parameters_;
diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl
index 54b98809b..4e1b591b2 100644
--- a/include/micm/solver/rosenbrock.inl
+++ b/include/micm/solver/rosenbrock.inl
@@ -433,8 +433,7 @@ namespace micm
 
   template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy, class LinearSolverPolicy>
   inline RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, LinearSolverPolicy>::RosenbrockSolver()
-      : system_(),
-        processes_(),
+      : processes_(),
         parameters_(RosenbrockSolverParameters::three_stage_rosenbrock_parameters()),
         state_parameters_(),
         process_set_(),
@@ -463,8 +462,7 @@ namespace micm
       const std::vector<Process>& processes,
       const RosenbrockSolverParameters& parameters,
       const std::function<LinearSolverPolicy(const SparseMatrixPolicy<double>, double)> create_linear_solver)
-      : system_(system),
-        processes_(processes),
+      : processes_(processes),
         parameters_(parameters),
         state_parameters_(),
         process_set_(),
@@ -474,7 +472,7 @@ namespace micm
     std::function<std::string(const std::vector<std::string>& variables, const std::size_t i)> state_reordering;
 
     std::size_t index = 0;
-    for (auto& name : system_.UniqueNames(state_reordering))
+    for (auto& name : system.UniqueNames())
       variable_map[name] = index++;
 
     // generate a state-vector reordering function to reduce fill-in in linear solver
@@ -483,7 +481,7 @@ namespace micm
       // get unsorted Jacobian non-zero elements
       auto unsorted_process_set = ProcessSet(processes, variable_map);
       auto unsorted_jac_elements = unsorted_process_set.NonZeroJacobianElements();
-      MatrixPolicy<int> unsorted_jac_non_zeros(system_.StateSize(), system_.StateSize(), 0);
+      MatrixPolicy<int> unsorted_jac_non_zeros(system.StateSize(), system.StateSize(), 0);
       for (auto& elem : unsorted_jac_elements)
         unsorted_jac_non_zeros[elem.first][elem.second] = 1;
       auto reorder_map = DiagonalMarkowitzReorder<MatrixPolicy>(unsorted_jac_non_zeros);
@@ -492,7 +490,7 @@ namespace micm
 
       variable_map.clear();
       std::size_t index = 0;
-      for (auto& name : system_.UniqueNames(state_reordering))
+      for (auto& name : system.UniqueNames(state_reordering))
         variable_map[name] = index++;
     }
     
@@ -508,7 +506,7 @@ namespace micm
     auto jacobian = build_jacobian<SparseMatrixPolicy>(
       process_set_.NonZeroJacobianElements(),
       parameters_.number_of_grid_cells_,
-      system_.StateSize()
+      system.StateSize()
     );
 
     std::vector<std::size_t> jacobian_diagonal_elements;
@@ -518,9 +516,10 @@ namespace micm
     state_parameters_ = {
       .number_of_grid_cells_ = parameters_.number_of_grid_cells_,
       .number_of_rate_constants_ = processes_.size(),
-      .variable_names_ = system_.UniqueNames(state_reordering),
+      .variable_names_ = system.UniqueNames(state_reordering),
       .custom_rate_parameter_labels_ = param_labels,
-      .jacobian_diagonal_elements_ = jacobian_diagonal_elements
+      .jacobian_diagonal_elements_ = jacobian_diagonal_elements,
+      .state_size_ = system.StateSize()
     };
 
     process_set_.SetJacobianFlatIds(jacobian);
@@ -535,7 +534,7 @@ namespace micm
     state.jacobian_ = build_jacobian<SparseMatrixPolicy>(
       process_set_.NonZeroJacobianElements(),
       state_parameters_.number_of_grid_cells_,
-      system_.StateSize()
+      state_parameters_.state_size_
     );
 
     auto lu = linear_solver_.GetLUMatrices(state.jacobian_, 1.0e-30);
diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp
index 46a0f3423..d6e112bd9 100644
--- a/include/micm/solver/state.hpp
+++ b/include/micm/solver/state.hpp
@@ -24,6 +24,7 @@ namespace micm
     std::vector<std::string> variable_names_{};
     std::vector<std::string> custom_rate_parameter_labels_{};
     std::vector<std::size_t> jacobian_diagonal_elements_;
+    size_t state_size_;
   };
 
   template<
diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl
index 5e8d7346f..27bd72b58 100644
--- a/include/micm/solver/state.inl
+++ b/include/micm/solver/state.inl
@@ -9,7 +9,8 @@ namespace micm
       : conditions_(),
         variables_(),
         custom_rate_parameters_(),
-        rate_constants_()
+        rate_constants_(),
+        jacobian_()
   {
   }
 
@@ -34,7 +35,8 @@ namespace micm
         rate_constants_(parameters.number_of_grid_cells_, parameters.number_of_rate_constants_, 0.0),
         variable_map_(),
         custom_rate_parameter_map_(),
-        variable_names_(parameters.variable_names_)
+        variable_names_(parameters.variable_names_),
+        jacobian_()
   {
     std::size_t index = 0;
     for (auto& name : variable_names_)
diff --git a/test/integration/e5.hpp b/test/integration/e5.hpp
index 3dbc37f17..efe968798 100644
--- a/test/integration/e5.hpp
+++ b/test/integration/e5.hpp
@@ -19,7 +19,6 @@ class E5 : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, Linea
      const micm::RosenbrockSolverParameters& parameters)
       : micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, LinearSolverPolicy>()
   {
-    this->system_ = system;
     this->processes_ = processes;
     this->parameters_ = parameters;
 
@@ -51,6 +50,7 @@ class E5 : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, Linea
       .variable_names_ = system.UniqueNames(state_reordering),
       .custom_rate_parameter_labels_ = param_labels,
       .jacobian_diagonal_elements_ = jacobian_diagonal_elements,
+      .state_size_ = system.StateSize()
     };
 
     this->linear_solver_ = LinearSolverPolicy(jacobian, 1.0e-30);
@@ -65,7 +65,7 @@ class E5 : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, Linea
     auto state = micm::State<MatrixPolicy, SparseMatrixPolicy>{ this->state_parameters_ };
 
     state.jacobian_ = micm::build_jacobian<SparseMatrixPolicy>(
-        nonzero_jacobian_elements_, this->state_parameters_.number_of_grid_cells_, this->system_.StateSize());
+        nonzero_jacobian_elements_, this->state_parameters_.number_of_grid_cells_, this->state_parameters_.state_size_);
 
     auto lu = this->linear_solver_.GetLUMatrices(state.jacobian_, 1.0e-30);
     auto lower_matrix = std::move(lu.first);
diff --git a/test/integration/hires.hpp b/test/integration/hires.hpp
index 5396fe4cf..4f7674bd8 100644
--- a/test/integration/hires.hpp
+++ b/test/integration/hires.hpp
@@ -20,7 +20,6 @@ class HIRES : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, Li
       const micm::RosenbrockSolverParameters& parameters)
       : micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, LinearSolverPolicy>()
   {
-    this->system_ = system;
     this->processes_ = processes;
     this->parameters_ = parameters;
 
@@ -52,6 +51,7 @@ class HIRES : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, Li
       .variable_names_ = system.UniqueNames(state_reordering),
       .custom_rate_parameter_labels_ = param_labels,
       .jacobian_diagonal_elements_ = jacobian_diagonal_elements,
+      .state_size_ = system.StateSize()
     };
 
     this->linear_solver_ = LinearSolverPolicy(jacobian, 1.0e-30);
@@ -66,7 +66,7 @@ class HIRES : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, Li
     auto state = micm::State<MatrixPolicy, SparseMatrixPolicy>{ this->state_parameters_ };
 
     state.jacobian_ = micm::build_jacobian<SparseMatrixPolicy>(
-        nonzero_jacobian_elements_, this->state_parameters_.number_of_grid_cells_, this->system_.StateSize());
+        nonzero_jacobian_elements_, this->state_parameters_.number_of_grid_cells_, this->state_parameters_.state_size_);
 
     auto lu = this->linear_solver_.GetLUMatrices(state.jacobian_, 1.0e-30);
     auto lower_matrix = std::move(lu.first);
diff --git a/test/integration/oregonator.hpp b/test/integration/oregonator.hpp
index 4f01511d7..e2477e0a8 100644
--- a/test/integration/oregonator.hpp
+++ b/test/integration/oregonator.hpp
@@ -20,7 +20,6 @@ class Oregonator : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolic
       const micm::RosenbrockSolverParameters& parameters)
       : micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, LinearSolverPolicy>()
   {
-    this->system_ = system;
     this->processes_ = processes;
     this->parameters_ = parameters;
 
@@ -52,6 +51,7 @@ class Oregonator : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolic
       .variable_names_ = system.UniqueNames(state_reordering),
       .custom_rate_parameter_labels_ = param_labels,
       .jacobian_diagonal_elements_ = jacobian_diagonal_elements,
+      .state_size_ = system.StateSize()
     };
 
     this->linear_solver_ = LinearSolverPolicy(jacobian, 1.0e-30);
@@ -66,7 +66,7 @@ class Oregonator : public micm::RosenbrockSolver<MatrixPolicy, SparseMatrixPolic
     auto state = micm::State<MatrixPolicy, SparseMatrixPolicy>{ this->state_parameters_ };
 
     state.jacobian_ = micm::build_jacobian<SparseMatrixPolicy>(
-        nonzero_jacobian_elements_, this->state_parameters_.number_of_grid_cells_, this->system_.StateSize());
+        nonzero_jacobian_elements_, this->state_parameters_.number_of_grid_cells_, this->state_parameters_.state_size_);
 
     auto lu = this->linear_solver_.GetLUMatrices(state.jacobian_, 1.0e-30);
     auto lower_matrix = std::move(lu.first);