diff --git a/mesons.jl b/mesons.jl
index 7c7fba9..243aa56 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -12,7 +12,7 @@ const g2_s = 110.349189097820 # dimensionless
 const g2_v = 187.694676506801 # dimensionless
 const g2_ρ = 192.927428365698 # dimensionless
 const g2_γ = 0.091701236 # dimensionless # equal to 4πα
-const κ = 3.260178893447 # MeV
+const κ_ss = 3.260178893447 # MeV
 const λ = -0.003551486718 # dimensionless # LambdaSS
 const ζ = 0.023499504053 # dimensionless # LambdaVV
 const Λv = 0.043376933644 # dimensionless # LambdaVR
@@ -26,15 +26,12 @@ greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * mi
 greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
-    m is the mass of the meson in MeV/c2,
-    r_max is the r-cutoff in fm."
-function solveKG(m, source, r_max)
-    Δr = r_max / (length(source) - 1)
-    rs = range(0, r_max; length=length(source))
-    int_measure = ħc .* Δr .* rs .^ 2
+    m is the mass of the meson in MeV/c2."
+function solveKG(m, source, s::system)
+    int_measure = ħc .* Δr(s) .* rs(s) .^ 2
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    greenMat = greensFunction.(rs, transpose(rs))
+    greenMat = greensFunction.(rs(s), transpose(rs(s)))
 
     return greenMat * (int_measure .* source)
 end
@@ -47,21 +44,21 @@ end
     r is the radius in fm,
     An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4))
+function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s::system, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4))
     (Φ0, W0, B0, A0) = initial_sol
-    (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
+    (src_Φ0, src_W0, src_B0, src_A0) = (zero_array(s) for _ in 1:4)
 
     # A0 doesn't need iterations
     @. src_A0 = -g2_γ * ρ_vp
-    A0 .= solveKG(m_γ, src_A0, r_max)
+    A0 .= solveKG(m_γ, src_A0, s)
 
     for _ in 1:iterations
-        @. src_Φ0 = g2_s * ((κ/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
+        @. 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, r_max)
-        W0 .= (solveKG(m_ω, src_W0, r_max) .+ W0) ./ 2 # to supress oscillation
-        B0 .= (solveKG(m_ρ, src_B0, r_max) .+ B0) ./ 2 # to supress oscillation
+        Φ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
     end
 
     return (Φ0, W0, B0, A0)
diff --git a/nucleons.jl b/nucleons.jl
index 32ce4a8..dc4fc2e 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -27,19 +27,19 @@ end
     divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
     shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
     the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
-    Δr = r_max / divs
-
+function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, s::system; shooting=true, normalize=true)
     if shooting
-        @assert divs % 2 == 0 "divs must be an even number when shooting=true"
-        prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2))
-        sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
+        @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
+        prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
+        sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
         wf_right = reverse(hcat(sol.u...); dims=2)
-        r_max = r_max / 2 # for the next step
+        next_r_max = s.r_max / 2 # for the next step
+    else
+        next_r_max = s.r_max
     end
     
-    prob = ODEProblem(dirac!, [0, 1], (0, r_max))
-    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
+    prob = ODEProblem(dirac!, [0, 1], (0, next_r_max))
+    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
     wf = hcat(sol.u...)
 
     if shooting # join two segments
@@ -52,7 +52,7 @@ function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, n
     end
 
     if normalize
-        norm = sum(wf .* wf) * Δr # integration by Reimann sum
+        norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
         wf = wf ./ sqrt(norm)
     end
 
@@ -62,30 +62,30 @@ end
 "Returns a function that solves the Dirac equation and returns g(r=r_max) where
     r_max is the outer boundary in fm,
     the other parameters are the same from dirac!(...)."
-function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=Tsit5())
-    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
-    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])[1, 1]
+function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, s::system; dtype=Float64, algo=Tsit5())
+    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
+    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1]
     return func
 end
 
 "Find all bound energies between E_min (=0) and E_max (=mass) where
     tol_digits determines the precision for root finding and the threshold for identifying duplicates,
     the other parameters are the same from dirac!(...)."
-function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
-    func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max)
+function findEs(κ, p, Φ0, W0, B0, A0, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
+    func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, s)
     Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
     return unique(E -> round(E; digits=tol_digits), Es)
 end
 
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
     the other parameters are defined above"
-function findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
+function findAllOrbitals(p, Φ0, W0, B0, A0, s::system, E_min=0, E_max=(p ? M_p : M_n))
     κs = Int[]
     Es = Float64[]
     # start from κ=-1 and go both up and down
     for direction in [-1, 1]
         for κ in direction * (1:100) # cutoff is 100
-            new_Es = findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min, E_max)
+            new_Es = findEs(κ, p, Φ0, W0, B0, A0, s, E_min, E_max)
             if isempty(new_Es); break; end
             append!(Es, new_Es)
             append!(κs, fill(κ, length(new_Es)))
@@ -94,20 +94,20 @@ function findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_
     return (κs, Es)
 end
 
-"For a given list of κ values with corresponding energies, attempt to fill N lowest lying orbitals and return occupancy numbers"
-function fillNucleons(N::Int, κs, Es)
+"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers"
+function fillNucleons(Z_or_N::Int, κs, Es)
     sort_i = sortperm(Es)
 
     occ = zeros(Int, length(κs))
 
     for i in sort_i
-        if N ≤ 0; break; end;
+        if Z_or_N ≤ 0; break; end;
         max_occ = 2 * j_κ(κs[i]) + 1
-        occ[i] = min(max_occ, N)
-        N -= occ[i]
+        occ[i] = min(max_occ, Z_or_N)
+        Z_or_N -= occ[i]
     end
 
-    N == 0 || @warn "All orbitals could not be filled"
+    Z_or_N == 0 || @warn "All orbitals could not be filled"
     return occ
 end
 
@@ -120,16 +120,16 @@ l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 "Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
     the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
     the other parameters are defined above"
-function calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, r_max, divs)
-    ρr2 = zeros(2, divs + 1) # ρ×r² values
+function calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, s::system)
+    ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
-        wf = solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
+        wf = solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, s; shooting=true, normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end
 
-    r2s = (collect ∘ range)(0, r_max, length=divs+1).^2
+    r2s = rs(s).^2
     ρ = ρr2 ./ transpose(r2s)
     ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
     
@@ -141,8 +141,8 @@ end
 
 "Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
     the other parameters are defined above"
-function solveNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=800, E_max=939)
-    κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min, E_max)
-    occs = fillNucleons(N, κs, Es)
-    return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, r_max, divs)
+function solveNucleonDensity(p, Φ0, W0, B0, A0, s::system, E_min=800, E_max=939)
+    κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, s, E_min, E_max)
+    occs = fillNucleons(Z_or_N(s, p), κs, Es)
+    return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, s)
 end
diff --git a/system.jl b/system.jl
index 4b425ef..8d9f82c 100644
--- a/system.jl
+++ b/system.jl
@@ -1,6 +1,4 @@
 using Interpolations, PolyLog, Plots
-include("nucleons.jl")
-include("mesons.jl")
 
 "Defines a nuclear system to be solved"
 struct system
@@ -11,13 +9,16 @@ struct system
 end
 
 "Get mass number of nucleus"
-A(s::system) = s.Z + s.N
+A(s::system)::Int = s.Z + s.N
 
 "Get r values in the mesh"
-rs(s::system) = range(0, s.r_max, length=s.divs+1)
+rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)
 
 "Get Δr value for the mesh"
-Δr(s::system) = s.r_max / s.divs
+Δr(s::system)::Float64 = s.r_max / s.divs
+
+"Get the number of protons or neutrons in the system"
+Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N
 
 "Create an empty array for the size of the mesh"
 zero_array(s::system) = zeros(1 + s.divs)
@@ -25,6 +26,9 @@ 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)))
 
+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)
@@ -46,7 +50,7 @@ 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.r_max, s.divs, isnan(E_total_previous) ? 50 : 5; initial_sol = (Φ0s, W0s, B0s, A0s))
+        @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)
@@ -54,14 +58,14 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
         A_interp = linear_interpolation(rs(s), A0s)
 
         # protons
-        @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, S_interp, V_interp, R_interp, A_interp, s.r_max)
+        @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.r_max, s.divs)
+        @time "Proton densities" (ρ_sp, ρ_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.r_max)
+        @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.r_max, s.divs)
+        @time "Neutron densities" (ρ_sn, ρ_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
diff --git a/test/Linear_nucleon_spectrum.jl b/test/Linear_nucleon_spectrum.jl
index 2f7b92f..e1b52fb 100644
--- a/test/Linear_nucleon_spectrum.jl
+++ b/test/Linear_nucleon_spectrum.jl
@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -16,14 +16,16 @@ R_interp = linear_interpolation(xs, Rs)
 A_interp = linear_interpolation(xs, As)
 
 p = false
-N = p ? 8 : 8
 r_max = maximum(xs)
+divs = length(xs) - 1
+s = system(8, 8, r_max, divs)
+
 E_min = 860
 E_max = 939
 
-κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(N, κs, Es)
+occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)
diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl
index e360d5f..9528951 100644
--- a/test/Pb208_meson_flds.jl
+++ b/test/Pb208_meson_flds.jl
@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../mesons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -23,8 +23,9 @@ xs = test_data[:, 1]
 
 r_max = maximum(xs)
 divs = length(xs) - 1
+s = system(0, 0, r_max, divs)
 
-(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 200) # many iterations needed without an initial solution
+(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, 200) # many iterations needed without an initial solution
 
 plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
 xlabel!("r (fm)")
\ No newline at end of file
diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl
index 481d55b..8588ee9 100644
--- a/test/Pb208_nucleon_basic.jl
+++ b/test/Pb208_nucleon_basic.jl
@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -18,15 +18,18 @@ A_interp = linear_interpolation(xs, As)
 κ = -1
 p = false
 r_max = maximum(xs)
+divs = length(xs) - 1
+s = system(0, 0, r_max, divs)
+
 E_min = 850
 E_max = 939
 
-boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
 
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
-func = boundaryValueFunc(κ, p, S_interp, V_interp, R_interp, A_interp, r_max)
+func = boundaryValueFunc(κ, p, S_interp, V_interp, R_interp, A_interp, s)
 Es = collect(E_min:0.5:E_max)
 boundaryVals = [func(E)^2 for E in Es]
 
diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl
index c99c23f..5c9c147 100644
--- a/test/Pb208_nucleon_dens.jl
+++ b/test/Pb208_nucleon_dens.jl
@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -18,21 +18,18 @@ A_interp = linear_interpolation(xs, As)
 N_p = 82
 N_n = 126
 r_max = maximum(xs)
-E_min = 860
-E_max = 939
-divs = 400
+divs = length(xs) - 1
+s = system(N_p, N_n, r_max, divs)
 
-rs = range(0, r_max, length=divs+1)
+(ρ_sp, ρ_vp) = solveNucleonDensity(true, S_interp, V_interp, R_interp, A_interp, s)
 
-(ρ_sp, ρ_vp) = solveNucleonDensity(N_p, true, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
+p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
-p_sp = plot(rs, ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(rs, ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+(ρ_sn, ρ_vn) = solveNucleonDensity(false, S_interp, V_interp, R_interp, A_interp, s)
 
-(ρ_sn, ρ_vn) = solveNucleonDensity(N_n, false, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
-
-p_sn = plot(rs, ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(rs, ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_vn = plot(rs(s), ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
 
 # benchmark data generated from Hartree.f
 # format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)
diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl
index 18992eb..5c8fa95 100644
--- a/test/Pb208_nucleon_spectrum.jl
+++ b/test/Pb208_nucleon_spectrum.jl
@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -16,14 +16,17 @@ R_interp = linear_interpolation(xs, Rs)
 A_interp = linear_interpolation(xs, As)
 
 p = true
-N = p ? 82 : 126
+N_p = 82
+N_n = 126
 r_max = maximum(xs)
-E_min = 860
+divs = length(xs) - 1
+s = system(N_p, N_n, r_max, divs)
+E_min = 800
 E_max = 939
 
-κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(N, κs, Es)
+occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)
diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl
index 0a13905..7724fee 100644
--- a/test/Pb208_nucleon_wf.jl
+++ b/test/Pb208_nucleon_wf.jl
@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -18,18 +18,18 @@ A_interp = linear_interpolation(xs, As)
 κ = -1
 p = true
 r_max = maximum(xs)
-E_min = 880
+divs = length(xs) - 1
+s = system(0, 0, r_max, divs)
+E_min = 800
 E_max = 939
 
-groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max) |> minimum
+groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max) |> minimum
 println("ground state E = $groundE")
 
-divs = 400
-wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
-rs = range(0, r_max, length=divs+1)
+wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, s)
 gs = wf[1, :]
 fs = wf[2, :]
 
-plot(rs, gs, label="g(r)")
-plot!(rs, fs, label="f(r)")
+plot(rs(s), gs, label="g(r)")
+plot!(rs(s), fs, label="f(r)")
 xlabel!("r (fm)")