diff --git a/src/harmonic.cpp b/src/harmonic.cpp
index 80ae318..8689505 100644
--- a/src/harmonic.cpp
+++ b/src/harmonic.cpp
@@ -2,6 +2,7 @@
 #include <igl/harmonic.h>
 #include <nanobind/nanobind.h>
 #include <nanobind/eigen/dense.h>
+#include <nanobind/eigen/sparse.h>
 
 namespace nb = nanobind;
 using namespace nb::literals;
@@ -23,6 +24,15 @@ namespace pyigl
     }
     return W;
   }
+  auto harmonic_integrated_from_laplacian_and_mass(
+    const Eigen::SparseMatrixN &L,
+    const Eigen::SparseMatrixN &M,
+    const int k)
+  {
+    Eigen::SparseMatrixN Q;
+    igl::harmonic(L, M, k, Q);
+    return Q;
+  }
 }
 
 // Bind the wrapper to the Python module
@@ -44,4 +54,17 @@ R"(Compute k-harmonic weight functions "coordinates".
 @param[in] bc #b by #W list of boundary values
 @param[in] k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
 @return W  #V by #W list of weights)");
-}
+  m.def(
+    "harmonic_integrated_from_laplacian_and_mass",
+    &pyigl::harmonic_integrated_from_laplacian_and_mass,
+    "L"_a, 
+    "M"_a, 
+    "k"_a, 
+R"(Build the discrete k-harmonic operator (computing integrated quantities). 
+That is, if the k-harmonic PDE is Q x = 0, then this minimizes x' Q x.
+
+@param[in] L  #V by #V discrete (integrated) Laplacian
+@param[in] M  #V by #V mass matrix
+@param[in] k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
+@return Q  #V by #V discrete (integrated) k-Laplacian)");
+}
\ No newline at end of file
diff --git a/tests/test_all.py b/tests/test_all.py
index d2067bc..7962135 100644
--- a/tests/test_all.py
+++ b/tests/test_all.py
@@ -157,6 +157,13 @@ def test_harmonic():
     W = igl.bbw(V,F,b,bc)
     W = igl.harmonic(V,F,b,bc,k=1)
     W = igl.harmonic(V,F,b,bc,k=2)
+
+def test_harmonic_integrated_from_laplacian_and_mass():
+    V, F = triangulated_square()
+    L = igl.cotmatrix(V, F)
+    M = igl.massmatrix(V, F, igl.MASSMATRIX_TYPE_VORONOI)
+    Q = igl.harmonic_integrated_from_laplacian_and_mass(L, M, k=1)
+    Q = igl.harmonic_integrated_from_laplacian_and_mass(L, M, k=2)
     
 def test_tets():
     V,F,T = single_tet()