diff --git a/nucleons.jl b/nucleons.jl
index 18d791d..1acdf96 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -28,18 +28,18 @@ end
     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)
-    saveat = r_max / divs
+    Δr = r_max / divs
 
     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=saveat)
+        sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
         wf_right = reverse(hcat(sol.u...); dims=2)
         r_max = r_max / 2 # for the next step
     end
     
     prob = ODEProblem(dirac!, [0, 1], (0, r_max))
-    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=saveat)
+    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
     wf = hcat(sol.u...)
 
     if shooting # join two segments
@@ -52,8 +52,8 @@ function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, n
     end
 
     if normalize
-        norm = sum(wf .* wf) * r_max / divs # integration by Reimann sum
-        wf = wf ./ sqrt(norm * 2) # WHY FACTOR OF 2?
+        norm = sum(wf .* wf) * Δr # integration by Reimann sum
+        wf = wf ./ sqrt(norm)
     end
 
     return wf