diff --git a/dirac.jl b/dirac.jl
index 590e0ad..a96307a 100644
--- a/dirac.jl
+++ b/dirac.jl
@@ -9,34 +9,34 @@ M_p = 938.2720813 # Proton mass in MeV/c2
 κ is the generalized angular momentum,
 M is the mass in MeV/c2,
 E in the energy in MeV,
-S(r) & V(r) are functions corresponding to scalar and vector potentials in MeV,
+Φ0(r) & W0(r) are functions corresponding to scalar and vector potentials in MeV,
 r is the radius in fm.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)."
-function dirac!(du, (g, f), (κ, M, E, S, V), r)
-    du[1] = -(κ/r) * g + (E + M - S(r) - V(r)) * f / ħc
-    du[2] =  (κ/r) * f - (E - M + S(r) - V(r)) * g / ħc
+function dirac!(du, (g, f), (κ, M, E, Φ0, W0), r)
+    du[1] = -(κ/r) * g + (E + M - Φ0(r) - W0(r)) * f / ħc
+    du[2] =  (κ/r) * f - (E - M + Φ0(r) - W0(r)) * g / ħc
 end
 
 "Solve the Dirac equation and return g(r=r_max) for given scalar and vector potentials where
 κ is the generalized angular momentum,
 M is the mass in MeV/c2,
 E in the energy in MeV,
-S(r) & V(r) are functions corresponding to scalar and vector potentials in MeV,
+Φ0(r) & W0(r) are functions corresponding to scalar and vector potentials in MeV,
 r_max is the outer boundary in fm,
 r_min (=r_max/1000) is inside boundary in fm which cannot be 0 due to the centrifugal term."
-function boundaryValue(κ, M, E, S, V, r_max, r_min=r_max/1000)
+function boundaryValue(κ, M, E, Φ0, W0, r_max, r_min=r_max/1000)
     prob = ODEProblem(dirac!, [0, 1], (r_min, r_max))
-    sol = solve(prob, RK4(), p=(κ, M, E, S, V))
+    sol = solve(prob, RK4(), p=(κ, M, E, Φ0, W0))
     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,
+Φ0(r) & W0(r) are functions corresponding to scalar and vector potentials in MeV,
 r_max is the outer boundary in fm,
 r_min (=r_max/1000) is inside boundary in fm 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)
+function findEs(κ, M, Φ0, W0, r_max, r_min=r_max/1000, E_min=0, E_max=M)
+    f(E) = boundaryValue(κ, M, E, Φ0, W0, r_max, r_min)
     return find_zeros(f, (E_min, E_max))
 end