diff --git a/system.jl b/system.jl
index 5b43948..02138b0 100644
--- a/system.jl
+++ b/system.jl
@@ -1,4 +1,4 @@
-using Interpolations, PolyLog
+using Interpolations, PolyLog, Plots
 include("nucleons.jl")
 include("mesons.jl")
 
@@ -26,7 +26,7 @@ zero_array(s::system) = zeros(1 + s.divs)
 Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a)))
 
 "Run the full program by self-consistent solution of nucleon and meson densities"
-function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4))
+function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4); monitor_print=true, monitor_plot=false)
     if isnothing(initial_dens)
         dens_guess = Woods_Saxon.(rs(s))
         ρ_sp = s.Z .* dens_guess
@@ -59,9 +59,11 @@ 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
+
         # total binding energy of nucleons
         E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n)
-        println("Total binding E per nucleon = $(E_total/A(s))")
+        monitor_print && println("Total binding E per nucleon = $(E_total/A(s))")
 
         # check convergence
         abs(E_total_previous - E_total) < 0.1 && break
diff --git a/test/Pb208.jl b/test/Pb208.jl
index b7da84e..e0b5d97 100644
--- a/test/Pb208.jl
+++ b/test/Pb208.jl
@@ -2,4 +2,4 @@ include("../system.jl")
 
 s = system(82, 126, 20.0, 400)
 
-solve_system(s)
+solve_system(s; monitor_plot=true)