diff --git a/src/cpp/divergence.cpp b/src/cpp/divergence.cpp
index f047d971..c1e97efb 100644
--- a/src/cpp/divergence.cpp
+++ b/src/cpp/divergence.cpp
@@ -27,8 +27,6 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       at(i, i - 1) = -1.0;
       at(i, i) = 1.0;
     }
-    // Weights
-    Q = { 1.0, 1.0, 1.0, 1.0, 1.0 };
     break;
   case 4:
     // A
@@ -50,9 +48,6 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       at(i, i) = 9.0 / 8.0;
       at(i, i + 1) = -1.0 / 24.0;
     }
-    Q  = { 2186.0 / 1943.0 , 2125.0 / 2828.0 , 1441.0 / 1240.0 , 648.0 / 673.0
-      , 349.0 / 350.0 , 648.0 / 673.0 , 1441.0 / 1240.0 , 2125.0 / 2828.0
-      , 2186.0 / 1943.0 };
     break;
   case 6:
     // A
@@ -94,11 +89,6 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       at(i, i + 1) = -25.0 / 384.0;
       at(i, i + 2) = 3.0 / 640.0;
     }
-    // Weights
-    Q  = { 2383.0 / 2005.0 , 929.0 / 2002.0 , 887.0 / 531.0 , 3124.0 / 5901.0
-      , 1706.0 / 1457.0 , 457.0 / 467.0 , 1057.0 / 1061.0 , 457.0 / 467.0
-      , 1706.0 / 1457.0 , 3124.0 / 5901.0 , 887.0 / 531.0 , 929.0 / 2002.0
-      , 2383.0 / 2005.0 };
     break;
   }
 
@@ -166,6 +156,3 @@ Divergence::Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
         Utils::spkron(A1, D1) + Utils::spkron(A2, D2) + Utils::spkron(A3, D3);
   }
 }
-
-// Returns weights
-vec Divergence::getQ() { return Q; }
diff --git a/src/cpp/divergence.h b/src/cpp/divergence.h
index 77e91538..9b599e99 100644
--- a/src/cpp/divergence.h
+++ b/src/cpp/divergence.h
@@ -61,15 +61,6 @@ class Divergence : public sp_mat {
    */  
   Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz);
   
-  /**
-   * @brief Returns the weights used in the Mimeitc Divergence Operators.
-   *
-   * @note for informational purposes only, can be used in constructing new operators.
-   */    
-  vec getQ();
-
-private:
-  vec Q;
 };
 
 #endif // DIVERGENCE_H
diff --git a/src/cpp/gradient.cpp b/src/cpp/gradient.cpp
index c2355a4c..6c38d67b 100644
--- a/src/cpp/gradient.cpp
+++ b/src/cpp/gradient.cpp
@@ -37,8 +37,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i) = -1.0;
       at(i, i + 1) = 1.0;
     }
-    // Weights
-    P = { 3.0 / 8.0 , 9.0 / 8.0 , 1.0 , 9.0 / 8.0 , 3.0 / 8.0 };
     break;
   case 4:
     // A
@@ -70,10 +68,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i + 1) = 9.0 / 8.0;
       at(i, i + 2) = -1.0 / 24.0;
     }
-    // Weights
-    P = { 1606.0 / 4535.0 , 941.0 / 766.0 , 1384.0 / 1541.0 , 1371.0 / 1346.0
-      , 701.0 / 700.0 , 1371.0 / 1346.0 , 1384.0 / 1541.0 , 941.0 / 766.0
-      , 1606.0 / 4535.0 };
     break;
   case 6:
     // A
@@ -129,12 +123,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i + 2) = -25.0 / 384.0;
       at(i, i + 3) = 3.0 / 640.0;
     }
-    // Weights
-    P = { 420249.0 / 1331069.0 , 2590978.0 / 1863105.0 , 882762.0 / 1402249.0
-      , 1677712.0 / 1359311.0 , 239985.0 / 261097.0 , 664189.0 / 657734.0
-      , 756049.0 / 754729.0 , 664189.0 / 657734.0 , 239985.0 / 261097.0
-      , 1677712.0 / 1359311.0 , 882762.0 / 1402249.0 , 2590978.0 / 1863105.0
-      , 420249.0 / 1331069.0 };
     break;
   case 8:
     // A
@@ -222,13 +210,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i + 3) = 49.0 / 5120.0;
       at(i, i + 4) = -5.0 / 7168.0;
     }
-    // Weights
-    P = { 267425.0 / 904736.0 , 2307435.0 / 1517812.0 , 847667.0 / 3066027.0
-      , 4050911.0 / 2301238.0 , 498943.0 / 1084999.0 , 211042.0 / 170117.0
-      , 2065895.0 / 2191686.0 , 1262499.0 / 1258052.0 , 1314891.0 / 1312727.0
-      , 1262499.0 / 1258052.0 , 2065895.0 / 2191686.0 , 211042.0 / 170117.0
-      , 498943.0 / 1084999.0 , 4050911.0 / 2301238.0 , 847667.0 / 3066027.0
-      , 2307435.0 / 1517812.0 , 267425.0 / 904736.0 };
     break;
   }
 
@@ -296,6 +277,3 @@ Gradient::Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
         Utils::spkron(A1, G1) + Utils::spkron(A2, G2) + Utils::spkron(A3, G3);
   }
 }
-
-// Returns weights
-vec Gradient::getP() { return P; }
diff --git a/src/cpp/gradient.h b/src/cpp/gradient.h
index 3b4f95a3..730a8aab 100644
--- a/src/cpp/gradient.h
+++ b/src/cpp/gradient.h
@@ -62,16 +62,6 @@ class Gradient : public sp_mat {
    */  
   Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz);
 
-
-  /**
-   * @brief Returns the weights used in the Mimeitc Gradient Operators.
-   *
-   * @note for informational purposes only, can be used in constructing new operators.
-   */  
-  vec getP();
-
-private:
-  vec P;
 };
 
 #endif // GRADIENT_H
diff --git a/src/cpp/mole.h b/src/cpp/mole.h
index 3fe0a3e1..fae7777c 100644
--- a/src/cpp/mole.h
+++ b/src/cpp/mole.h
@@ -24,5 +24,7 @@
 #include "operators.h"
 #include "robinbc.h"
 #include "utils.h"
+#include "weightsQ.h"
+#include "weightsP.h"
 
 #endif // MOLE_H
diff --git a/src/cpp/utils.cpp b/src/cpp/utils.cpp
index b6f02d6d..4642cc86 100644
--- a/src/cpp/utils.cpp
+++ b/src/cpp/utils.cpp
@@ -18,8 +18,9 @@
 #include "utils.h"
 #include <cassert>
 
-#ifdef EIGEN
+//#ifdef EIGEN
 #include <eigen3/Eigen/SparseLU>
+#include <eigen3/Eigen/SparseQR>
 
 vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) {
   Eigen::SparseMatrix<Real> eigen_A(A.n_rows, A.n_cols);
@@ -47,7 +48,34 @@ vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) {
 
   return vec(eigen_x.data(), eigen_x.size());
 }
-#endif
+
+vec Utils::spsolve_eigenQR(const sp_mat &A, const vec &b) {
+  Eigen::SparseMatrix<Real> eigen_A(A.n_rows, A.n_cols);
+  std::vector<Eigen::Triplet<Real>> triplets;
+  Eigen::SparseQR<Eigen::SparseMatrix<Real>, Eigen::COLAMDOrdering<int>> solver;
+
+  Eigen::VectorXd eigen_x(A.n_rows);
+  triplets.reserve(5 * A.n_rows);
+
+  auto it = A.begin();
+  while (it != A.end()) {
+    triplets.push_back(Eigen::Triplet<Real>(it.row(), it.col(), *it));
+    ++it;
+  }
+
+  eigen_A.setFromTriplets(triplets.begin(), triplets.end());
+  triplets.clear();
+
+  auto b_ = conv_to<std::vector<Real>>::from(b);
+  Eigen::Map<Eigen::VectorXd> eigen_b(b_.data(), b_.size());
+
+  solver.analyzePattern(eigen_A);
+  solver.factorize(eigen_A);
+  eigen_x = solver.solve(eigen_b);
+
+  return vec(eigen_x.data(), eigen_x.size());
+}
+//#endif
 
 // Basic implementation of Kronecker product
 /*
diff --git a/src/cpp/utils.h b/src/cpp/utils.h
index dc860a85..e3fb27ac 100644
--- a/src/cpp/utils.h
+++ b/src/cpp/utils.h
@@ -68,6 +68,8 @@ class Utils {
   */
   static vec spsolve_eigen(const sp_mat &A, const vec &b);
 
+  static vec spsolve_eigenQR(const sp_mat &A, const vec &b);
+
   /**
   * @brief An analog to the MATLAB/Octave 2D meshgrid operation
   *
diff --git a/src/cpp/weightsP.cpp b/src/cpp/weightsP.cpp
new file mode 100644
index 00000000..05b1c3c7
--- /dev/null
+++ b/src/cpp/weightsP.cpp
@@ -0,0 +1,32 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+
+/*
+ * @file weights.cpp
+ * 
+ * @brief Mimetic Differences Q Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#include "utils.h"
+#include "divergence.h"
+#include "gradient.h"
+#include "weightsP.h"
+
+WeightsP::WeightsP(u16 k, u32 m, Real dx)
+{
+  Gradient G(k, m, dx);
+
+  arma::vec b(m+2);
+  b.at(0) = -1.0;
+  b.at(m+1) = 1.0;
+
+  sp_mat Gtranspose = G.t();
+  *this = Utils::spsolve_eigenQR(Gtranspose,b);
+
+}
\ No newline at end of file
diff --git a/src/cpp/weightsP.h b/src/cpp/weightsP.h
new file mode 100644
index 00000000..ea490978
--- /dev/null
+++ b/src/cpp/weightsP.h
@@ -0,0 +1,42 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+ 
+ /*
+ * @file weights.h
+ * 
+ * @brief Mimetic Differences Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#ifndef WEIGHTSP_H
+#define WEIGHTSP_H
+
+#include "utils.h"
+#include <cassert>
+
+/**
+ * @brief Generate P Weights
+ *
+ */
+class WeightsP : public sp_vec {
+
+public:
+  using sp_vec::operator=;
+
+  /**
+   * @brief P Weights Constructor
+   *
+   * @param k Order of accuracy
+   * @param m Number of cells
+   * @param dx Spacing between cells
+   */  
+  WeightsP(u16 k, u32 m, Real dx);
+
+};
+
+#endif // WEIGHTSP_H
diff --git a/src/cpp/weightsQ.cpp b/src/cpp/weightsQ.cpp
new file mode 100644
index 00000000..23bef9dc
--- /dev/null
+++ b/src/cpp/weightsQ.cpp
@@ -0,0 +1,43 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+
+/*
+ * @file weights.cpp
+ * 
+ * @brief Mimetic Differences Q Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#include "utils.h"
+#include "divergence.h"
+#include "gradient.h"
+#include "weightsQ.h"
+
+WeightsQ::WeightsQ(u16 k, u32 m, Real dx): sp_vec(m + 1) 
+{
+  Divergence D(k, m, dx);
+  D.shed_row(0);
+  D.shed_row(m);
+
+  vec b(m+1);
+  b.at(0) = -1.0;
+  b.at(m) = 1.0;
+
+  sp_mat Dtranspose = D.t();
+  sp_vec Q_prime = Utils::spsolve_eigenQR(Dtranspose, b);
+  sp_vec Q(m+2);
+  Q.at(0) = 1.0;
+  Q.at(m+1) = 1.0;
+  for (int i = 1; i < m+1; ++i) {
+    Q.at(i) = Q_prime.at(i-1);
+  }
+
+  *this = Q;
+
+}
+
diff --git a/src/cpp/weightsQ.h b/src/cpp/weightsQ.h
new file mode 100644
index 00000000..597f92b4
--- /dev/null
+++ b/src/cpp/weightsQ.h
@@ -0,0 +1,42 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+ 
+ /*
+ * @file weights.h
+ * 
+ * @brief Mimetic Differences Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#ifndef WEIGHTSQ_H
+#define WEIGHTSQ_H
+
+#include "utils.h"
+#include <cassert>
+
+/**
+ * @brief Generate Q Weights
+ *
+ */
+class WeightsQ : public sp_vec {
+
+public:
+  using sp_vec::operator=;
+
+  /**
+   * @brief Q Weights Constructor
+   *
+   * @param k Order of accuracy
+   * @param m Number of cells
+   * @param dx Spacing between cells
+   */  
+  WeightsQ(u16 k, u32 m, Real dx);
+
+};
+
+#endif // WEIGHTSQ_H
diff --git a/tests/cpp/test6.cpp b/tests/cpp/test6.cpp
new file mode 100644
index 00000000..208bd19c
--- /dev/null
+++ b/tests/cpp/test6.cpp
@@ -0,0 +1,66 @@
+#include "mole.h"
+#include <gtest/gtest.h>
+
+arma::vec Q_baseline[3];
+arma::vec P_baseline[3];
+constexpr Real tolerance = 1.0e-08;
+
+void run_Qtest(int k) {
+    Real dx = 1.0;
+    int m = k*2 + 1;
+    int j = k/2 - 1;
+    WeightsQ Q(k,m,dx);
+    Real total_error = 0.0;
+
+    for (int i = 0; i < m+2; ++i) {
+        total_error += abs(Q[i] - Q_baseline[j][i]);
+    }
+    ASSERT_LE(total_error, tolerance);
+}
+
+void run_Ptest(int k) {
+    Real dx = 1.0;
+    int m = k*2 + 1;
+    int j = k/2 - 1;
+    Real total_error = 0.0;
+
+    WeightsP P(k,m,dx);
+    for (int i = 0; i < m; ++i) {
+        total_error += abs(P[i] - P_baseline[j][i]);
+    }
+    ASSERT_LE(total_error, tolerance);
+
+}
+
+TEST(WeightTests, Accuracy) {
+
+        // Baseline weights - Generated from MatLab/Octave
+        Q_baseline[0] = { 
+            1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000
+        };
+        Q_baseline[1] = {
+            1.000000000,1.125064293,0.751414447,1.162097097,0.962852897,0.997142531,
+            0.962852897,1.162097097,0.751414447,1.125064293,1.000000000
+        };
+        Q_baseline[2] = {
+            1.000000000,1.188528786,0.464036031,1.670433241,0.529401810,1.170898848,0.978586058,0.996230453,
+            0.978586058,1.170898848,0.529401810,1.670433241,0.464036031,1.188528786,1.000000000
+        };
+
+        P_baseline[0] = { 
+            0.375000000,1.125000000,1.000000000,1.000000000,1.125000000,0.375000000
+        };
+        P_baseline[1] = {
+            0.354134518,1.228459404,0.898117099,1.018547095,1.000741884,
+            1.000741884,1.018547095,0.898117099,1.228459404,0.354134518
+        };
+        P_baseline[2] = {
+            0.315722926,1.390677390,0.629532984,1.234237262,0.919144328,1.009804636,1.000880474,
+            1.000880474,1.009804636,0.919144328,1.234237262,0.629532984,1.390677390,0.315722926
+        };
+      
+        for (int k : {2,4,6}) {
+            run_Qtest(k);
+            run_Ptest(k);
+        }
+}