Fixed the unstability

mesons.jl

@@ -41,7 +41,7 @@ end
     r is the radius in fm,
     the inital solutions are read from s and the final solutions are saved in-place.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function solveMesonFields!(s::system, iterations=10)
+function solveMesonFields!(s::system, iterations=15, oscillation_control_parameter=0.3)
     (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
     (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
 
@@ -55,9 +55,15 @@ function solveMesonFields!(s::system, iterations=10)
         @. src_Φ0 = g2_s * ((κ_ss/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
         @. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (2B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
         @. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2
+
         Φ0 .= solveKG(m_s, src_Φ0, s)
-        W0 .= (solveKG(m_ω, src_W0, s) .+ W0) ./ 2 # to supress oscillation
-        B0 .= (solveKG(m_ρ, src_B0, s) .+ B0) ./ 2 # to supress oscillation
+
+        # W0 and B0 keep a fraction of their previous solutions to suppress oscillations
+        W0 .*= (1 - oscillation_control_parameter)
+        B0 .*= (1 - oscillation_control_parameter)
+
+        W0 .+= solveKG(m_ω, src_W0, s) .* oscillation_control_parameter
+        B0 .+= solveKG(m_ρ, src_B0, s) .* oscillation_control_parameter
     end
 
     return (Φ0, W0, B0, A0)

system.jl

@@ -73,7 +73,7 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
 
     while true
         # mesons
-        @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 5)
+        @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
 
         # protons
         @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s)