diff --git a/mesons.jl b/mesons.jl
index e2c2312..85c6f0c 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -42,10 +42,14 @@ 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)
-    (Φ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)
+    # A0 doesn't need iterations
+    src_A0 = ρ_vp ./ ħc
+    A0 = solveKG(0, src_A0, r_max)
 
-    for _ in iterations
+    (Φ0, W0, B0) = (zeros(1 + divs) for _ in 1:3)
+    (src_Φ0, src_W0, src_B0) = (zeros(1 + divs) for _ in 1:3)
+
+    for _ in 1:iterations
         @. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn)) / ħc
         Φ0 .= solveKG(m_s / sqrt(ħc), src_Φ0, r_max)
 
@@ -53,10 +57,7 @@ function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3)
         W0 .= solveKG(m_ω / sqrt(ħc), src_W0, r_max)
 
         @. src_B0 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2) / ħc
-        W0 .= solveKG(m_ρ / sqrt(ħc), src_B0, r_max)
-
-        @. src_A0 = ρ_vp / ħc
-        A0 .= solveKG(0, src_A0, r_max) # this doesn't need iterations
+        B0 .= solveKG(m_ρ / sqrt(ħc), src_B0, r_max)
     end
 
     return (Φ0, W0, B0, A0)