Root finding implemented

Project.toml

@@ -1,3 +1,4 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"

dirac.jl

@@ -1,4 +1,4 @@
-using DifferentialEquations
+using DifferentialEquations, Roots
 
 ħc = 197.327 # ħc in MeVfm
 M_n = 939.5654133 # Neutron mass in MeV/c2
@@ -29,3 +29,14 @@ function boundaryValue(κ, M, E, S, V, r_max, r_min=r_max/1000)
     sol = solve(prob, RK4(), p=(κ, M, E, S, V))
     return sol(r_max)[1]
 end
+
+"Find all bound energies between E_min (=0) and E_max (=M) where
+κ is the generalized angular momentum,
+M is the mass in MeV/c2,
+S(r) & V(r) are functions corresponding to scalar and vector potentials in MeV,
+r_max is the outer boundary,
+r_min (=r_max/1000) is inside boundary which cannot be 0 due to the centrifugal term."
+function findEs(κ, M, S, V, r_max, r_min=r_max/1000, E_min=0, E_max=M)
+    f(E) = boundaryValue(κ, M, E, S, V, r_max, r_min)
+    return find_zeros(f, (E_min, E_max))
+end

test/Pb208.jl

@@ -15,8 +15,18 @@ V_interp = linear_interpolation(xs, Vs)
 R_interp = linear_interpolation(xs, Rs)
 A_interp = linear_interpolation(xs, As)
 
-Es = collect(840:0.5:940)
-boundaryVals = [boundaryValue(-1, M_n, E, S_interp, V_interp, maximum(xs))^2 for E in Es]
+κ = -1
+M = M_n
+r_max = maximum(xs)
+E_min = M - 100
+E_max = M
+
+boundEs = findEs(κ, M_n, S_interp, V_interp, r_max, r_max/1000, E_min, E_max)
+println("bound E = $boundEs")
+
+Es = collect(E_min:0.5:E_max)
+boundaryVals = [boundaryValue(κ, M_n, E, S_interp, V_interp, r_max)^2 for E in Es]
 
 plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
+vline!(boundEs, label="bound E")
 xlabel!("E (MeV)")