diff --git a/system.jl b/system.jl
index 02138b0..181892e 100644
--- a/system.jl
+++ b/system.jl
@@ -39,6 +39,10 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
 
     (Φ0s, W0s, B0s, A0s) = initial_flds
 
+    if monitor_plot
+        p = plot(legends=false, layout=(2, 2), title=["ρ_sp" "ρ_vp" "ρ_sn" "ρ_vn"])
+    end
+
     E_total_previous = NaN
 
     while true
@@ -59,7 +63,13 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
         occs_n = fillNucleons(s.N, κs_n, Es_n)
         (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
 
-        monitor_plot && plot(rs(s), hcat(ρ_sp, ρ_vp, ρ_sn, ρ_vn), layout=(2, 2), title=["ρ_sp" "ρ_vp" "ρ_sn" "ρ_vn"]) |> display
+        if monitor_plot
+            for s in p.series_list
+                s.plotattributes[:linecolor] = :gray
+            end
+            plot!(p, rs(s), hcat(ρ_sp, ρ_vp, ρ_sn, ρ_vn), linecolor=:red)
+            display(p)
+        end
 
         # total binding energy of nucleons
         E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n)