diff --git a/mesons.jl b/mesons.jl
index 1481acd..b669710 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -25,7 +25,10 @@ greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
     m is the mass of the meson in MeV/c2."
 function solveKG(m, source, s::system)
-    int_measure = ħc .* Δr(s) .* rs(s) .^ 2
+
+    @assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
+    simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
+    int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
     greenMat = greensFunction.(rs(s), transpose(rs(s)))