diff --git a/main.cpp b/main.cpp
index 0c9d875..ae8ed68 100644
--- a/main.cpp
+++ b/main.cpp
@@ -51,6 +51,9 @@ struct FEM {
   std::unique_ptr<mfem::GridFunction> H_gf;
   std::unique_ptr<mfem::GridFunction> P_gf;
 
+  mfem::Vector com;
+  mfem::DenseMatrix Q;
+
   [[nodiscard]] bool okay() const;
 };
 
@@ -214,7 +217,11 @@ void write_output(const Envelope& env, std::ostream& out, const Args& args);
 template <typename T>
 std::map<std::string, T> invert_pair_map(const std::map<T, std::pair<const char *, const char *>>& forward_map);
 
+mfem::Vector get_com(const FEM& fem, const mfem::GridFunction &rho);
 
+mfem::DenseMatrix compute_quadrupole_moment_tensor(const FEM& fem, const mfem::GridFunction& rho, const mfem::Vector& com);
+
+double l2_multipole_potential(const FEM& fem, const mfem::Vector& x);
 
 
 /*****************************
@@ -889,3 +896,97 @@ void write_output(const Envelope& env, std::ostream& out, const Args& args) {
   }
 }
 
+mfem::Vector get_com(const FEM& fem, const mfem::GridFunction &rho) {
+  const int dim = fem.mesh_ptr->Dimension();
+  mfem::Vector com(dim);
+  com = 0.0;
+  double total_mass = 0.0;
+
+  for (int i = 0; i < fem.H1_fes->GetNE(); ++i) {
+    mfem::ElementTransformation *trans = fem.H1_fes->GetElementTransformation(i);
+    const mfem::IntegrationRule &ir = mfem::IntRules.Get(trans->GetGeometryType(), fem.H1_fes->GetOrder(0) + trans->OrderW());
+
+    for (int j = 0; j < ir.GetNPoints(); ++j) {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(j);
+      trans->SetIntPoint(&ip);
+
+      double weight = trans->Weight() * ip.weight;
+      double rho_val = rho.GetValue(i, ip);
+
+      mfem::Vector phys_point(dim);
+      trans->Transform(ip, phys_point);
+
+      const double mass_term = rho_val * weight;
+      total_mass += mass_term;
+
+      for (int d = 0; d < dim; ++d) {
+        com(d) += phys_point(d) * mass_term;
+      }
+    }
+  }
+  com /= total_mass;
+  return com;
+}
+
+mfem::DenseMatrix compute_quadrupole_moment_tensor(const FEM& fem, const mfem::GridFunction& rho, const mfem::Vector& com) {
+  const int dim = fem.mesh_ptr->Dimension();
+  mfem::DenseMatrix Q(dim, dim);
+  Q = 0.0;
+
+  for (int i = 0; i < fem.H1_fes->GetNE(); ++i) {
+    mfem::ElementTransformation *trans = fem.mesh_ptr->GetElementTransformation(i);
+    const mfem::IntegrationRule &ir = mfem::IntRules.Get(trans->GetGeometryType(), 2 * fem.H1_fes->GetOrder(0) + trans->OrderW());
+    for (int j = 0; j < ir.GetNPoints(); ++j) {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(j);
+      trans->SetIntPoint(&ip);
+
+      const double weight = trans->Weight() * ip.weight;
+      const double rho_val = rho.GetValue(i, ip);
+
+      mfem::Vector phys_point(dim);
+      trans->Transform(ip, phys_point);
+
+      mfem::Vector x_prime(dim);
+      double r_sq = 0.0;
+
+      for (int d = 0; d < dim; ++d) {
+        x_prime(d) = phys_point(d) - com(d);
+        r_sq += x_prime(d) * x_prime(d);
+      }
+
+      for (int m = 0; m < dim; ++m) {
+        for (int n = 0; n < dim; ++n) {
+          const double delta = (m == n) ? 1.0 : 0.0;
+          const double contrib = 3.0 * x_prime(m) * x_prime(n) - delta * r_sq;
+          Q(m, n) += rho_val * contrib * weight;
+        }
+      }
+    }
+  }
+  return Q;
+}
+
+double l2_multipole_potential(const FEM &fem, const mfem::Vector &x) {
+  const double r = x.Norml2();
+  if (r < 1e-12) return 0.0;
+
+  const int dim = fem.mesh_ptr->Dimension();
+
+  mfem::Vector n(x);
+  n /= r;
+
+  double l2_mult_factor = 0.0;
+  for (int i = 0; i < dim; ++i) {
+    for (int j = 0; j < dim; ++j) {
+      l2_mult_factor += fem.Q(i, j) * n(i) * n(j);
+    }
+  }
+
+  const double l2_contrib = (G / (2.0 * std::pow(r, 3))) * l2_mult_factor;
+
+  const double total_mass = get_current_mass(fem, *fem.rho_gf);
+  const double l0_contrib = -G * total_mass / r;
+
+  // l1 contribution is zero for a system centered on its COM
+  return l0_contrib + l2_contrib;
+}