diff --git a/system.jl b/system.jl
index 8d9f82c..4490727 100644
--- a/system.jl
+++ b/system.jl
@@ -1,11 +1,28 @@
 using Interpolations, PolyLog, Plots
 
-"Defines a nuclear system to be solved"
-struct system
+"Defines a nuclear system containing relevant parameters and meson/nucleon densities"
+mutable struct system
     Z::Int
     N::Int
+
     r_max::Float64
     divs::Int
+
+    Φ0::Vector{Float64}
+    W0::Vector{Float64}
+    B0::Vector{Float64}
+    A0::Vector{Float64}
+
+    ρ_sp::Vector{Float64}
+    ρ_vp::Vector{Float64}
+    ρ_sn::Vector{Float64}
+    ρ_vn::Vector{Float64}
+
+    "Dummy struct to define the mesh"
+    system(r_max, divs) = new(0, 0, r_max, divs, [Float64[] for _ in 1:8]...)
+
+    "Initialize an unsolved system"
+    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...)
 end
 
 "Get mass number of nucleus"
@@ -26,23 +43,28 @@ zero_array(s::system) = zeros(1 + s.divs)
 "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)))
 
+"Create linear interpolations for the meson fields"
+function fields_interp(s::system)
+    S_interp = linear_interpolation(rs(s), s.Φ0)
+    V_interp = linear_interpolation(rs(s), s.W0)
+    R_interp = linear_interpolation(rs(s), s.B0)
+    A_interp = linear_interpolation(rs(s), s.A0)
+    return (S_interp, V_interp, R_interp, A_interp)
+end
+
 include("nucleons.jl")
 include("mesons.jl")
 
 "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); monitor_print=true, monitor_plot=false)
-    if isnothing(initial_dens)
+function solve_system(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false)
+    if reinitialize_densities
         dens_guess = Woods_Saxon.(rs(s))
-        ρ_sp = s.Z .* dens_guess
-        ρ_vp = s.Z .* dens_guess
-        ρ_sn = s.N .* dens_guess
-        ρ_vn = s.N .* dens_guess
-    else
-        (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = initial_dens
+        @. s.ρ_sp = s.Z * dens_guess
+        @. s.ρ_vp = s.Z * dens_guess
+        @. s.ρ_sn = s.N * dens_guess
+        @. s.ρ_vn = s.N * dens_guess
     end
 
-    (Φ0s, W0s, B0s, A0s) = initial_flds
-
     if monitor_plot
         p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
     end
@@ -50,28 +72,25 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
     E_total_previous = NaN
 
     while true
-        @time "Meson fields" (Φ0s, W0s, B0s, A0s) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, isnan(E_total_previous) ? 50 : 5; initial_sol = (Φ0s, W0s, B0s, A0s))
-
-        S_interp = linear_interpolation(rs(s), Φ0s)
-        V_interp = linear_interpolation(rs(s), W0s)
-        R_interp = linear_interpolation(rs(s), B0s)
-        A_interp = linear_interpolation(rs(s), A0s)
+        # mesons
+        @time "Meson fields" solveMesonWfs(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s, isnan(E_total_previous) ? 50 : 5; initial_sol = (s.Φ0, s.W0, s.B0, s.A0))
+        (S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
 
         # protons
         @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, S_interp, V_interp, R_interp, A_interp, s)
         occs_p = fillNucleons(s.Z, κs_p, Es_p)
-        @time "Proton densities" (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s)
+        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s)
 
         # neutrons
         @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, S_interp, V_interp, R_interp, A_interp, s)
         occs_n = fillNucleons(s.N, κs_n, Es_n)
-        @time "Neutron densities" (ρ_sn, ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s)
+        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s)
 
         if monitor_plot
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
-            plot!(p, rs(s), hcat(ρ_sp, ρ_vp, ρ_sn, ρ_vn, Φ0s, W0s, B0s, A0s), linecolor=:red)
+            plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
             display(p)
         end
 
diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl
index 9528951..ffbf48f 100644
--- a/test/Pb208_meson_flds.jl
+++ b/test/Pb208_meson_flds.jl
@@ -23,7 +23,7 @@ xs = test_data[:, 1]
 
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(0, 0, r_max, divs)
+s = system(r_max, divs)
 
 (Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, 200) # many iterations needed without an initial solution
 
diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl
index 8588ee9..6349a39 100644
--- a/test/Pb208_nucleon_basic.jl
+++ b/test/Pb208_nucleon_basic.jl
@@ -1,4 +1,4 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
 include("../system.jl")
 
 # test data generated from Hartree.f
@@ -10,16 +10,18 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 κ = -1
 p = false
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(0, 0, r_max, divs)
+s = system(r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
+(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
 
 E_min = 850
 E_max = 939
diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl
index 5c9c147..6ffb7a9 100644
--- a/test/Pb208_nucleon_dens.jl
+++ b/test/Pb208_nucleon_dens.jl
@@ -1,4 +1,4 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
 include("../system.jl")
 
 # test data generated from Hartree.f
@@ -10,17 +10,19 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
 s = system(N_p, N_n, r_max, divs)
 
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
+(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
+
 (ρ_sp, ρ_vp) = solveNucleonDensity(true, S_interp, V_interp, R_interp, A_interp, s)
 
 p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl
index 5c8fa95..ed20529 100644
--- a/test/Pb208_nucleon_spectrum.jl
+++ b/test/Pb208_nucleon_spectrum.jl
@@ -10,17 +10,20 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 p = true
 N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
 s = system(N_p, N_n, r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
+(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
+
 E_min = 800
 E_max = 939
 
diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl
index 7724fee..59543f4 100644
--- a/test/Pb208_nucleon_wf.jl
+++ b/test/Pb208_nucleon_wf.jl
@@ -10,16 +10,19 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 κ = -1
 p = true
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(0, 0, r_max, divs)
+s = system(r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
+(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
+
 E_min = 800
 E_max = 939