diff --git a/mesons.jl b/mesons.jl
index 217be65..0295aa1 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -3,19 +3,19 @@ using DifferentialEquations
 const ħc = 197.327 # MeVfm
 
 # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
-# Values taken from Hartree.f (FSUGarnet) and
-const m_ρ = 763.0 # MeV/c2
-const m_ω = 782.5 # MeV/c2
+# Values taken from Hartree.f (FSUGarnet)
 const m_s = 496.939473213388 # MeV/c2
-const m_γ = 0.000001000 # MeV/c2 # not defined in paper
-const g2_s = 110.349189097820 # units?
-const g2_v = 187.694676506801 # units?
-const g2_ρ = 192.927428365698 # units?
+const m_ω = 782.5 # MeV/c2
+const m_ρ = 763.0 # MeV/c2
+const m_γ = 0.000001000 # MeV/c2 # should be 0?
+const g2_s = 110.349189097820 # dimensionless
+const g2_v = 187.694676506801 # dimensionless
+const g2_ρ = 192.927428365698 # dimensionless
 const g2_γ = 0.091701236 # dimensionless # equal to 4πα
-const κ = 3.260178893447 # units?
-const λ = -0.003551486718 # units? # LambdaSS
-const ζ = 0.023499504053 # units? # LambdaVV
-const Λv = 0.043376933644 # units? # LambdaVR
+const κ = 3.260178893447
+const λ = -0.003551486718 # LambdaSS
+const ζ = 0.023499504053 # LambdaVV
+const Λv = 0.043376933644 # LambdaVR
 
 const r_reg = 1E-9 # fm # regulator for Green's functions
 
@@ -47,21 +47,21 @@ end
     r is the radius in fm.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
 function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3)
-    # A0 doesn't need iterations
-    src_A0 = -g2_γ .* ρ_vp .* ħc
-    A0 = solveKG(0, src_A0, r_max)
+    (Φ0, W0, B0, A0) = (zeros(1 + divs) for _ in 1:4)
+    (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
 
-    (Φ0, W0, B0) = (zeros(1 + divs) for _ in 1:3)
-    (src_Φ0, src_W0, src_B0) = (zeros(1 + divs) for _ in 1:3)
+    # A0 doesn't need iterations
+    @. src_A0 = -g2_γ * ρ_vp * ħc
+    A0 .= solveKG(m_γ, src_A0, r_max)
 
     for _ in 1:iterations
-        @. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn)) * ħc
+        @. src_Φ0 = g2_s * ((κ/2) * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn) * ħc
         Φ0 .= solveKG(m_s, src_Φ0, r_max)
 
-        @. src_W0 = g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp + ρ_vn)) * ħc
+        @. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2 * Λv * (B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn) * ħc
         W0 .= solveKG(m_ω, src_W0, r_max)
 
-        @. src_B0 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2) * ħc
+        @. src_B0 = g2_ρ * (2 * Λv * (W0/ħc)^2 * (B0/ħc)) - (g2_ρ / 2) * (ρ_vp - ρ_vn) * ħc
         B0 .= solveKG(m_ρ, src_B0, r_max)
     end