diff --git a/poet/include/IPhreeqcPOET.hpp b/poet/include/IPhreeqcPOET.hpp
index 0b66bf90..1a550242 100644
--- a/poet/include/IPhreeqcPOET.hpp
+++ b/poet/include/IPhreeqcPOET.hpp
@@ -14,6 +14,7 @@
 #include <map>
 #include <string>
 #include <sys/types.h>
+#include <utility>
 #include <vector>
 
 enum { POET_SOL = 0, POET_EXCH, POET_KIN, POET_EQUIL, POET_SURF };
@@ -33,42 +34,51 @@ public:
     this->LoadDatabase(database.c_str());
     this->RunFile(input_script.c_str());
 
-    this->parseInitValues();
+    this->parseInit();
   }
 
-  std::vector<std::string> getInitNames() const {
-    std::vector<std::string> names;
+  IPhreeqcPOET(const std::string &database, const std::string &input_script,
+               const std::vector<std::string> &solutionInitVector,
+               std::uint32_t n_cells)
+      : IPhreeqc(), n_cells(n_cells), solutionInitVector(solutionInitVector) {
+    this->LoadDatabase(database.c_str());
+    this->RunFile(input_script.c_str());
 
-    for (auto &sol : this->initial_names) {
-      names.insert(names.end(), sol.begin(), sol.end());
+    if (n_cells > 1) {
+      std::string copy_string =
+          "COPY cell 1 2-" + std::to_string(n_cells) + "\n";
     }
-
-    return names;
   }
 
-  std::vector<std::vector<double>> getInitValues() const {
-    return this->initial_concentrations;
-  }
+  struct PhreeqcMat {
+    std::vector<std::string> names;
+    std::vector<int> ids;
+    std::vector<std::vector<double>> values;
+  };
 
-  std::vector<int> getSolutionIds() const { return this->solution_ids; }
+  PhreeqcMat getPhreeqcMat();
 
+  PhreeqcMat getPhreeqcMat(const std::vector<int> &ids);
 
   std::map<int, std::string> raw_dumps() {
     std::map<int, std::string> dumps;
 
     this->SetDumpStringOn(true);
 
-    for (const auto &sol_id : this->solution_ids) {
+    for (const auto &[sol_id, unused] : this->raw_initials) {
       std::string call_string = "DUMP\n -cells " + std::to_string(sol_id);
       this->RunString(call_string.c_str());
       dumps[sol_id] = this->GetDumpString();
     }
 
+    this->SetDumpStringOn(false);
+
     return dumps;
   }
 
+  using essential_names = std::array<std::vector<std::string>, 5>;
   using ModulesArray = std::array<std::uint32_t, 5>;
-  
+
   ModulesArray getModuleSizes() const {
     ModulesArray module_sizes;
     for (std::uint8_t i = 0; i < 5; i++) {
@@ -79,13 +89,13 @@ public:
   }
 
 private:
-  using essential_names = std::array<std::vector<std::string>, 5>;
-
+  // required only for simulation
   essential_names initial_names;
-  std::vector<std::vector<double>> initial_concentrations;
-  std::vector<int> solution_ids;
+  std::uint32_t n_cells;
+  std::vector<std::string> solutionInitVector;
 
   void parseInitValues();
+  void parseInit();
 
   essential_names dump_essential_names(std::size_t cell_number);
 
@@ -110,13 +120,15 @@ private:
     return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n);
   }
 
-  void resolveSolutionUseKW(
-      const std::vector<IPhreeqc::SolutionMapping> &unresolved,
-      std::map<int, std::pair<essential_names, std::vector<double>>>
-          &mapped_values);
+  using RawMap = std::map<int, std::pair<essential_names, std::vector<double>>>;
 
-  void valuesFromModule(const std::string &module_name, int cell_number,
-                        essential_names &names, std::vector<double> &values);
+  essential_names union_raws(const RawMap &raws);
+
+  std::vector<std::vector<double>>
+  conc_from_essentials(const RawMap &raws, const essential_names &names);
+
+  // void valuesFromModule(const std::string &module_name, int cell_number,
+  //                       essential_names &names, std::vector<double> &values);
 
   std::string subExchangeName(std::string name) {
     for (const auto &species : this->PhreeqcPtr->Get_species_list()) {
@@ -128,4 +140,6 @@ private:
     }
     return name;
   };
+
+  RawMap raw_initials;
 };
\ No newline at end of file
diff --git a/poet/src/GetSet.cpp b/poet/src/GetSet.cpp
new file mode 100644
index 00000000..8a7c4aff
--- /dev/null
+++ b/poet/src/GetSet.cpp
@@ -0,0 +1,119 @@
+#include <IPhreeqcPOET.hpp>
+
+IPhreeqcPOET::essential_names
+IPhreeqcPOET::dump_essential_names(std::size_t cell_number) {
+
+  essential_names eNames;
+
+  // Solutions
+  if (this->Get_solution(cell_number) != NULL) {
+    std::vector<std::string> &eSolNames = eNames[POET_SOL];
+    this->Get_solution(cell_number)->dump_essential_names(eSolNames);
+  }
+
+  // Exchange
+  if (this->Get_exchange(cell_number) != NULL) {
+    std::vector<std::string> &eExchNames = eNames[POET_EXCH];
+    this->Get_exchange(cell_number)->dump_essential_names(eExchNames);
+  }
+
+  // Kinetics
+  if (this->Get_kinetic(cell_number) != NULL) {
+    std::vector<std::string> &eKinNames = eNames[POET_KIN];
+    this->Get_kinetic(cell_number)->dump_essential_names(eKinNames);
+  }
+
+  // PPassemblage
+  if (this->Get_equilibrium(cell_number) != NULL) {
+    std::vector<std::string> &eEquNames = eNames[POET_EQUIL];
+    this->Get_equilibrium(cell_number)->dump_essential_names(eEquNames);
+  }
+
+  // Surface
+  if (this->Get_surface(cell_number) != NULL) {
+    std::vector<std::string> &eSurfNames = eNames[POET_SURF];
+    this->Get_surface(cell_number)->dump_essential_names(eSurfNames);
+  }
+
+  return eNames;
+}
+
+std::vector<double>
+IPhreeqcPOET::get_essential_values(std::size_t cell_number,
+                                   const std::vector<std::string> &order) {
+  std::vector<double> essentials;
+
+  // Solutions
+  if (this->Get_solution(cell_number) != NULL) {
+    std::vector<double> sol_values;
+
+    this->Get_solution(cell_number)->get_essential_values(sol_values, order);
+    essentials.insert(essentials.end(), sol_values.begin(), sol_values.end());
+  }
+
+  // Exchange
+  if (this->Get_exchange(cell_number) != NULL) {
+    std::vector<double> exch_values;
+
+    this->Get_exchange(cell_number)->get_essential_values(exch_values);
+    essentials.insert(essentials.end(), exch_values.begin(), exch_values.end());
+  }
+
+  // Kinetics
+  if (this->Get_kinetic(cell_number) != NULL) {
+    std::vector<double> kin_values;
+
+    this->Get_kinetic(cell_number)->get_essential_values(kin_values);
+    essentials.insert(essentials.end(), kin_values.begin(), kin_values.end());
+  }
+
+  // PPassemblage
+  if (this->Get_equilibrium(cell_number) != NULL) {
+    std::vector<double> equ_values;
+
+    this->Get_equilibrium(cell_number)->get_essential_values(equ_values);
+    essentials.insert(essentials.end(), equ_values.begin(), equ_values.end());
+  }
+
+  // Surface
+  if (this->Get_surface(cell_number) != NULL) {
+    std::vector<double> surf_values;
+
+    this->Get_surface(cell_number)->get_essential_values(surf_values);
+    essentials.insert(essentials.end(), surf_values.begin(), surf_values.end());
+  }
+
+  return essentials;
+}
+
+void IPhreeqcPOET::set_essential_values(std::size_t cell_number,
+                                        const std::vector<std::string> &order,
+                                        std::vector<double> &values) {
+
+  auto dump_it = values.begin();
+
+  // Solutions
+  if (this->Get_solution(cell_number) != NULL) {
+    this->Get_solution(cell_number)->set_essential_values(dump_it, order);
+  }
+
+  // Exchange
+  if (this->Get_exchange(cell_number) != NULL) {
+    this->Get_exchange(cell_number)->set_essential_values(dump_it);
+  }
+
+  // Kinetics
+  if (this->Get_kinetic(cell_number) != NULL) {
+    this->Get_kinetic(cell_number)->set_essential_values(dump_it);
+  }
+
+  // PPassemblage
+  if (this->Get_equilibrium(cell_number) != NULL) {
+    this->Get_equilibrium(cell_number)->set_essential_values(dump_it);
+  }
+
+  // Surface
+  if (this->Get_surface(cell_number) != NULL) {
+    this->Get_surface(cell_number)->set_essential_values(dump_it);
+  }
+}
\ No newline at end of file
diff --git a/poet/src/IPhreeqcPOET.cpp b/poet/src/IPhreeqcPOET.cpp
index a64a6be6..099c0301 100644
--- a/poet/src/IPhreeqcPOET.cpp
+++ b/poet/src/IPhreeqcPOET.cpp
@@ -38,57 +38,58 @@ createConcVector(const std::vector<std::string> &conc_names,
   return conc_vec;
 }
 
-void IPhreeqcPOET::valuesFromModule(const std::string &module_name,
-                                    int cell_number, essential_names &names,
-                                    std::vector<double> &values) {
-  std::size_t dest_module_i = 0;
-  std::vector<double> to_insert;
-  if (module_name == "exchange") { // 1
-    this->Get_exchange(cell_number)->dump_essential_names(names[POET_EXCH]);
-    this->Get_exchange(cell_number)->get_essential_values(to_insert);
-    dest_module_i = 1;
-  } else if (module_name == "kinetics") { // 2
-    this->Get_kinetic(cell_number)->dump_essential_names(names[POET_KIN]);
-    this->Get_kinetic(cell_number)->get_essential_values(to_insert);
-    dest_module_i = 2;
-  } else if (module_name == "surface") { // 4
-    this->Get_surface(cell_number)->dump_essential_names(names[POET_SURF]);
-    this->Get_surface(cell_number)->get_essential_values(to_insert);
-    dest_module_i = 4;
-  }
+// void IPhreeqcPOET::valuesFromModule(const std::string &module_name,
+//                                     int cell_number, essential_names &names,
+//                                     std::vector<double> &values) {
+//   std::size_t dest_module_i = 0;
+//   std::vector<double> to_insert;
+//   if (module_name == "exchange") { // 1
+//     this->Get_exchange(cell_number)->dump_essential_names(names[POET_EXCH]);
+//     this->Get_exchange(cell_number)->get_essential_values(to_insert);
+//     dest_module_i = 1;
+//   } else if (module_name == "kinetics") { // 2
+//     this->Get_kinetic(cell_number)->dump_essential_names(names[POET_KIN]);
+//     this->Get_kinetic(cell_number)->get_essential_values(to_insert);
+//     dest_module_i = 2;
+//   } else if (module_name == "surface") { // 4
+//     this->Get_surface(cell_number)->dump_essential_names(names[POET_SURF]);
+//     this->Get_surface(cell_number)->get_essential_values(to_insert);
+//     dest_module_i = 4;
+//   }
 
-  std::size_t offset = 0;
-  for (std::size_t i = 0; i < dest_module_i; i++) {
-    offset += names[i].size();
-  }
+//   std::size_t offset = 0;
+//   for (std::size_t i = 0; i < dest_module_i; i++) {
+//     offset += names[i].size();
+//   }
 
-  values.insert(values.begin() + offset, to_insert.begin(), to_insert.end());
-}
+//   values.insert(values.begin() + offset, to_insert.begin(), to_insert.end());
+// }
 
-void IPhreeqcPOET::resolveSolutionUseKW(
-    const std::vector<IPhreeqc::SolutionMapping> &unresolved,
-    std::map<int, std::pair<essential_names, std::vector<double>>>
-        &mapped_values) {
+// void IPhreeqcPOET::resolveSolutionUseKW(
+//     const std::vector<IPhreeqc::SolutionMapping> &unresolved,
+//     std::map<int, std::pair<essential_names, std::vector<double>>>
+//         &mapped_values) {
 
-  for (const auto &input : unresolved) {
-    if (mapped_values.find(input.module_n) != mapped_values.end()) {
-      continue;
-    }
+//   for (const auto &input : unresolved) {
+//     if (mapped_values.find(input.module_n) != mapped_values.end()) {
+//       continue;
+//     }
 
-    essential_names new_conc_names;
-    new_conc_names[0] = mapped_values[input.sol_n].first[0];
+//     essential_names new_conc_names;
+//     new_conc_names[0] = mapped_values[input.sol_n].first[0];
 
-    const auto &curr_sol_vec = mapped_values[input.sol_n].second;
-    std::vector<double> new_conc_values(
-        curr_sol_vec.begin(), curr_sol_vec.begin() + new_conc_names[0].size());
+//     const auto &curr_sol_vec = mapped_values[input.sol_n].second;
+//     std::vector<double> new_conc_values(
+//         curr_sol_vec.begin(), curr_sol_vec.begin() +
+//         new_conc_names[0].size());
 
-    valuesFromModule(input.module_name, input.module_n, new_conc_names,
-                     new_conc_values);
+//     valuesFromModule(input.module_name, input.module_n, new_conc_names,
+//                      new_conc_values);
 
-    mapped_values[input.module_n] =
-        std::make_pair(new_conc_names, new_conc_values);
-  }
-}
+//     mapped_values[input.module_n] =
+//         std::make_pair(new_conc_names, new_conc_values);
+//   }
+// }
 
 void IPhreeqcPOET::parseInitValues() {
 
@@ -106,8 +107,8 @@ void IPhreeqcPOET::parseInitValues() {
         curr_conc_names, this->get_essential_values(id, curr_conc_names[0]));
   }
 
-  const auto unresolved_modules = this->getSolutionMapping();
-  resolveSolutionUseKW(unresolved_modules, init_values);
+  // const auto unresolved_modules = this->getSolutionMapping();
+  // resolveSolutionUseKW(unresolved_modules, init_values);
 
   std::vector<int> ids_to_erase;
 
@@ -182,122 +183,4 @@ void IPhreeqcPOET::parseInitValues() {
   //         solution_only, init_values[val.sol_n].second);
   //   }
   // }
-}
-
-IPhreeqcPOET::essential_names
-IPhreeqcPOET::dump_essential_names(std::size_t cell_number) {
-
-  essential_names eNames;
-
-  // Solutions
-  if (this->Get_solution(cell_number) != NULL) {
-    std::vector<std::string> &eSolNames = eNames[POET_SOL];
-    this->Get_solution(cell_number)->dump_essential_names(eSolNames);
-  }
-
-  // Exchange
-  if (this->Get_exchange(cell_number) != NULL) {
-    std::vector<std::string> &eExchNames = eNames[POET_EXCH];
-    this->Get_exchange(cell_number)->dump_essential_names(eExchNames);
-  }
-
-  // Kinetics
-  if (this->Get_kinetic(cell_number) != NULL) {
-    std::vector<std::string> &eKinNames = eNames[POET_KIN];
-    this->Get_kinetic(cell_number)->dump_essential_names(eKinNames);
-  }
-
-  // PPassemblage
-  if (this->Get_equilibrium(cell_number) != NULL) {
-    std::vector<std::string> &eEquNames = eNames[POET_EQUIL];
-    this->Get_equilibrium(cell_number)->dump_essential_names(eEquNames);
-  }
-
-  // Surface
-  if (this->Get_surface(cell_number) != NULL) {
-    std::vector<std::string> &eSurfNames = eNames[POET_SURF];
-    this->Get_surface(cell_number)->dump_essential_names(eSurfNames);
-  }
-
-  return eNames;
-}
-
-std::vector<double>
-IPhreeqcPOET::get_essential_values(std::size_t cell_number,
-                                   const std::vector<std::string> &order) {
-  std::vector<double> essentials;
-
-  // Solutions
-  if (this->Get_solution(cell_number) != NULL) {
-    std::vector<double> sol_values;
-
-    this->Get_solution(cell_number)->get_essential_values(sol_values, order);
-    essentials.insert(essentials.end(), sol_values.begin(), sol_values.end());
-  }
-
-  // Exchange
-  if (this->Get_exchange(cell_number) != NULL) {
-    std::vector<double> exch_values;
-
-    this->Get_exchange(cell_number)->get_essential_values(exch_values);
-    essentials.insert(essentials.end(), exch_values.begin(), exch_values.end());
-  }
-
-  // Kinetics
-  if (this->Get_kinetic(cell_number) != NULL) {
-    std::vector<double> kin_values;
-
-    this->Get_kinetic(cell_number)->get_essential_values(kin_values);
-    essentials.insert(essentials.end(), kin_values.begin(), kin_values.end());
-  }
-
-  // PPassemblage
-  if (this->Get_equilibrium(cell_number) != NULL) {
-    std::vector<double> equ_values;
-
-    this->Get_equilibrium(cell_number)->get_essential_values(equ_values);
-    essentials.insert(essentials.end(), equ_values.begin(), equ_values.end());
-  }
-
-  // Surface
-  if (this->Get_surface(cell_number) != NULL) {
-    std::vector<double> surf_values;
-
-    this->Get_surface(cell_number)->get_essential_values(surf_values);
-    essentials.insert(essentials.end(), surf_values.begin(), surf_values.end());
-  }
-
-  return essentials;
-}
-
-void IPhreeqcPOET::set_essential_values(std::size_t cell_number,
-                                        const std::vector<std::string> &order,
-                                        std::vector<double> &values) {
-
-  auto dump_it = values.begin();
-
-  // Solutions
-  if (this->Get_solution(cell_number) != NULL) {
-    this->Get_solution(cell_number)->set_essential_values(dump_it, order);
-  }
-
-  // Exchange
-  if (this->Get_exchange(cell_number) != NULL) {
-    this->Get_exchange(cell_number)->set_essential_values(dump_it);
-  }
-
-  // Kinetics
-  if (this->Get_kinetic(cell_number) != NULL) {
-    this->Get_kinetic(cell_number)->set_essential_values(dump_it);
-  }
-
-  // PPassemblage
-  if (this->Get_equilibrium(cell_number) != NULL) {
-    this->Get_equilibrium(cell_number)->set_essential_values(dump_it);
-  }
-
-  // Surface
-  if (this->Get_surface(cell_number) != NULL) {
-    this->Get_surface(cell_number)->set_essential_values(dump_it);
-  }
 }
\ No newline at end of file
diff --git a/poet/src/Init.cpp b/poet/src/Init.cpp
new file mode 100644
index 00000000..6ff4e80f
--- /dev/null
+++ b/poet/src/Init.cpp
@@ -0,0 +1,110 @@
+#include <IPhreeqcPOET.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <tuple>
+#include <vector>
+
+static inline std::vector<std::string>
+unionStringVectors(const std::vector<std::string> &a,
+                   const std::vector<std::string> &b) {
+  std::vector<std::string> result;
+
+  std::set_union(a.begin(), a.end(), b.begin(), b.end(),
+                 std::back_inserter(result));
+
+  return result;
+}
+
+static std::vector<double>
+createConcVector(const std::vector<std::string> &conc_names,
+                 const std::vector<double>::const_iterator &conc_it,
+                 const std::vector<std::string> &new_names) {
+  std::vector<double> conc_vec;
+
+  for (const auto &i : new_names) {
+    auto it = std::find(conc_names.begin(), conc_names.end(), i);
+
+    auto conc_index = std::distance(conc_names.begin(), it);
+
+    if (it != conc_names.end()) {
+      conc_vec.push_back(*(conc_it + conc_index));
+    } else {
+      conc_vec.push_back(std::numeric_limits<double>::quiet_NaN());
+    }
+  }
+
+  return conc_vec;
+}
+
+void IPhreeqcPOET::parseInit() {
+
+  std::map<int, std::pair<essential_names, std::vector<double>>> init_values;
+
+  for (const auto &[id, val] : this->PhreeqcPtr->Get_Rxn_solution_map()) {
+    // A key less than zero indicates an internal solution
+    if (id < 0) {
+      continue;
+    }
+
+    const essential_names curr_conc_names = this->dump_essential_names(id);
+
+    auto pair = std::make_pair(
+        curr_conc_names, this->get_essential_values(id, curr_conc_names[0]));
+    raw_initials[id] = pair;
+  }
+}
+
+IPhreeqcPOET::essential_names IPhreeqcPOET::union_raws(const RawMap &raws) {
+  IPhreeqcPOET::essential_names names;
+  for (const auto &[sol_id, val] : raws) {
+    for (std::size_t i = 0; i < names.size(); i++) {
+      names[i] = unionStringVectors(names[i], val.first[i]);
+    }
+
+    for (auto &specie : names[POET_EXCH]) {
+      specie = this->subExchangeName(specie);
+    }
+  }
+
+  return names;
+}
+
+std::vector<std::vector<double>>
+IPhreeqcPOET::conc_from_essentials(const RawMap &raws,
+                                   const essential_names &names) {
+  std::vector<std::vector<double>> values;
+
+  for (const auto &[key, val] : raws) {
+    std::vector<double> combined_conc_vec;
+
+    for (std::size_t i = 0, offset = 0; i < names.size(); i++) {
+      std::vector<double> union_vec =
+          createConcVector(val.first[i], val.second.begin() + offset, names[i]);
+
+      combined_conc_vec.insert(combined_conc_vec.end(), union_vec.begin(),
+                               union_vec.end());
+
+      offset += val.first[i].size();
+    }
+
+    values.push_back(combined_conc_vec);
+  }
+}
+
+IPhreeqcPOET::PhreeqcMat IPhreeqcPOET::getPhreeqcMat() {
+  std::vector<std::vector<double>> values;
+
+  const IPhreeqcPOET::essential_names names =
+      this->union_raws(this->raw_initials);
+
+  const std::vector<std::vector<double>> conc_values =
+      this->conc_from_essentials(this->raw_initials, names);
+
+
+  std::vector<int> ids(this->raw_initials.size());
+
+  std::transform(this->raw_initials.begin(), this->raw_initials.end(),
+                 ids.begin(), [](const auto &pair) { return pair.first; });
+
+}