diff --git a/nucleons.jl b/nucleons.jl
index 7fbc194..3c3f81d 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -51,6 +51,9 @@ function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64)
     return [g, f]
 end
 
+"Initial boundary condition for u(r)=[g(r), f(r)] at r=0"
+init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]?
+
 "Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
     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,
@@ -68,7 +71,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
         next_r_max = s.r_max
     end
     
-    prob = ODEProblem(dirac!, [0, 1], (0, next_r_max))
+    prob = ODEProblem(dirac!, init_BC(), (0, next_r_max))
     sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
     wf = hcat(sol.u...)
 
@@ -97,7 +100,7 @@ end
     where is r is in fm."
 function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
-    prob_left  = ODEProblem(dirac!, [0.0, 1.0], (0, r))
+    prob_left  = ODEProblem(dirac!, init_BC(), (0, r))
     function func(E)
         prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r))
         u_left  = solve(prob_left,  algo, p=(κ, E, f1, f2), saveat=[r])