diff --git a/system.jl b/system.jl
index 181892e..af14994 100644
--- a/system.jl
+++ b/system.jl
@@ -61,7 +61,7 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
         # neutrons
         κs_n, Es_n = findAllOrbitals(false, S_interp, V_interp, R_interp, A_interp, s.r_max)
         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)
+        (ρ_sn, ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
 
         if monitor_plot
             for s in p.series_list