diff --git a/NuclearRMF.jl b/NuclearRMF.jl
index 5d3f331..c71d400 100644
--- a/NuclearRMF.jl
+++ b/NuclearRMF.jl
@@ -11,8 +11,19 @@ end
 "Normalized Woods-Saxon form used for constructing an initial solution"
 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)))
 
+"Conditionally toggle @time if enable_time is true"
+macro conditional_time(label, expr)
+    quote
+        if $(esc(:enable_time))
+            @time $label $(esc(expr))
+        else
+            $(esc(expr))
+        end
+    end
+end
+
 "Run the full program by self-consistent solution of nucleon and meson densities"
-function solve_system!(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false)
+function solve_system!(s::system; reinitialize_densities=true, print_E=true, print_time=false, live_plots=false)
     if reinitialize_densities
         dens_guess = Woods_Saxon.(rs(s))
         @. s.ρ_sp = s.Z * dens_guess
@@ -21,18 +32,19 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
         @. s.ρ_vn = s.N * dens_guess
     end
 
-    if monitor_plot
+    if live_plots
         p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
     end
 
     previous_E_per_A = NaN
 
     while true
-        @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
-        @time "Proton densities" solveNucleonDensity!(true, s)
-        @time "Neutron densities" solveNucleonDensity!(false, s)
+        enable_time = print_time
+        @conditional_time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
+        @conditional_time "Proton densities" solveNucleonDensity!(true, s)
+        @conditional_time "Neutron densities" solveNucleonDensity!(false, s)
 
-        if monitor_plot
+        if live_plots
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
@@ -41,7 +53,7 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
         end
 
         E_per_A = total_E(s) / A(s)
-        monitor_print && println("Total binding E per nucleon = $E_per_A")
+        print_E && println("Total binding E per nucleon = $E_per_A")
 
         # check convergence
         abs(previous_E_per_A - E_per_A) < 0.0001 && break
diff --git a/test/Pb208.jl b/test/Pb208.jl
index cae50b9..38fe28a 100644
--- a/test/Pb208.jl
+++ b/test/Pb208.jl
@@ -7,4 +7,4 @@ divs = 400
 
 s = system(Z, N, r_max, divs)
 
-solve_system!(s; monitor_plot=true)
+solve_system!(s; live_plots=false)