diff --git a/mesons.jl b/mesons.jl
index b669710..96f90de 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -71,14 +71,14 @@ end
 
 "Calculate the total energy associated with meson fields"
 function meson_E(s::system)
-    int = 0.0
-    for (r, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(rs(s), s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+    E_densities = map(zip(s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)) do (Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn)
         E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
         E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
         E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
         E_γ = -(1/2) * (A0/ħc) * ρ_vp
         E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
-        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2
+        E_σ + E_ω + E_ρ + E_γ + E_ωρ
     end
-    return 4π * int * Δr(s) * ħc
+    
+    return simpsons_integrate(E_densities .* rs(s).^2, Δr(s); coefficient=4π * ħc)
 end