diff --git a/mesons.jl b/mesons.jl
index 75bf59a..7ac95af 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -20,19 +20,23 @@ const Λv = 0.043376933644 # units? # LambdaVR
 const r_reg = 1E-9 # fm # regulator for Green's functions
 
 "Green's function for Klein-Gordon equation in natural units"
-greensFunctionKG(m, r, rp) = -(1 / m) * rp / (r + r_reg) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
+greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
 
 "Green's function for Poisson's equation in natural units"
-greensFunctionP(r, rp) = -rp^2 / (max(r, rp) + r_reg)
+greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) for a given a source function (as an array) where
     m is the mass of the meson in MeV/c2,
     r_max is the r-cutoff in fm."
 function solveKG(m, source, r_max)
+    Δr = r_max / (length(source) - 1)
+    rs = range(0, r_max; length=length(source))
+    int_measure = Δr .* rs .^ 2
+
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    mesh = range(0, r_max; length=length(source))
-    greenMat = greensFunction.(mesh, transpose(mesh))
-    return greenMat * source
+    greenMat = greensFunction.(rs, transpose(rs))
+
+    return greenMat * (int_measure .* source)
 end
 
 "Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where