diff --git a/common.jl b/common.jl
index ccb93b5..c0b7503 100644
--- a/common.jl
+++ b/common.jl
@@ -1,2 +1,17 @@
 const ħc = 197.33 # MeVfm
 const r_reg = 1E-8 # fm # regulator for R
+
+"Integrate a uniformly discretized function f using Simpson's rule where h is the step size and coefficient is an optional scaling factor."
+function simpsons_integrate(f::AbstractVector{Float64}, h::Float64; coefficient::Float64 = 1.0)
+    @assert length(f) % 2 == 1 "Number of mesh divisions must be even for Simpson's rule"
+    s = sum(enumerate(f)) do (i, fi)
+        if i == 1 || i == length(f)
+            return fi
+        elseif i % 2 == 0
+            return 4fi
+        else
+            return 2fi
+        end
+    end
+    return (h/3) * coefficient * s
+end
diff --git a/nucleons.jl b/nucleons.jl
index 393637c..63b3d83 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -89,7 +89,9 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
-        wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum
+        g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
+        f2_int = simpsons_integrate(wf[2, :] .^ 2, Δr(s))
+        wf ./= sqrt(g2_int + f2_int)
     end
 
     return wf