diff --git a/nucleons.jl b/nucleons.jl
index 095c212..564534a 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -177,8 +177,10 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
 
     r2s = s.r_mesh.r .^ 2
     ρ = ρr2 ./ transpose(r2s)
-    ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
-    
+    if !all(isfinite.(ρ[:, 1]))
+        ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
+    end
+
     ρ_s = ρ[1, :] - ρ[2, :] # g^2 - f^2
     ρ_v = ρ[1, :] + ρ[2, :] # g^2 + f^2