diff --git a/NuclearRMF.jl b/NuclearRMF.jl
index 9877aad..602e1ef 100644
--- a/NuclearRMF.jl
+++ b/NuclearRMF.jl
@@ -11,7 +11,7 @@ Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli
 "Run the full program by self-consistent solution of nucleon and meson densities"
 function solve_system!(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false)
     if reinitialize_densities
-        dens_guess = Woods_Saxon.(rs(s))
+        dens_guess = Woods_Saxon.(s.r_mesh.r)
         @. s.ρ_sp = s.Z * dens_guess
         @. s.ρ_vp = s.Z * dens_guess
         @. s.ρ_sn = s.N * dens_guess
@@ -33,7 +33,7 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
-            plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
+            plot!(p, s.r_mesh.r, hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
             display(p)
         end
 
diff --git a/mesons.jl b/mesons.jl
index 1481acd..7829421 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -25,16 +25,15 @@ greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
     m is the mass of the meson in MeV/c2."
 function solveKG(m, source, s::system)
-    int_measure = ħc .* Δr(s) .* rs(s) .^ 2
+    int_measure = ħc .* s.r_mesh.w .* (s.r_mesh.r .^ 2)
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    greenMat = greensFunction.(rs(s), transpose(rs(s)))
+    greenMat = greensFunction.(s.r_mesh.r, transpose(s.r_mesh.r))
 
     return greenMat * (int_measure .* source)
 end
 
 "Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where
-    divs is the number of mesh divisions so the arrays are of length (1+divs),
     r is the radius in fm,
     the inital solutions are read from s and the final solutions are saved in-place.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
@@ -69,13 +68,13 @@ end
 "Calculate the total energy associated with meson fields"
 function meson_E(s::system)
     int = 0.0
-    for (r, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(rs(s), s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+    for (r, w, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(s.r_mesh.r, s.r_mesh.w, s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
         E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
         E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
         E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
         E_γ = -(1/2) * (A0/ħc) * ρ_vp
         E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
-        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2
+        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2 * w
     end
-    return 4π * int * Δr(s) * ħc
+    return 4π * int * ħc
 end
diff --git a/nucleons.jl b/nucleons.jl
index e57e0f5..d42a52f 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -33,44 +33,45 @@ function optimized_dirac_potentials(p::Bool, s::system)
     @. f1s =  M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
     @. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
 
-    f1 = linear_interpolation(rs(s), f1s)
-    f2 = linear_interpolation(rs(s), f2s)
+    f1 = linear_interpolation(s.r_mesh.r, f1s)
+    f2 = linear_interpolation(s.r_mesh.r, f2s)
 
     return (f1, f2)
 end
 
 "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,
+    the solution would be returned as a 2×mesh_size matrix,
     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::Bool, E, s::system; shooting=true, normalize=true, algo=Tsit5())
+function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Tsit5())
     (f1, f2) = optimized_dirac_potentials(p, s)
 
-    if shooting
-        @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
-        prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
-        sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
-        wf_right = reverse(hcat(sol.u...); dims=2)
-        next_r_max = s.r_max / 2 # for the next step
-    else
-        next_r_max = s.r_max
-    end
-    
-    prob = ODEProblem(dirac!, [0, 1], (0, next_r_max))
-    sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
-    wf = hcat(sol.u...)
+    # partitioning
+    mid_idx = length(s.r_mesh.r) ÷ 2
+    r_mid = s.r_mesh.r[mid_idx]
+    left_r = s.r_mesh.r[1:mid_idx]
+    right_r = s.r_mesh.r[mid_idx:end]
 
-    if shooting # join two segments
-        rescale_factor_g = wf[1, end] / wf_right[1, 1]
-        rescale_factor_f = wf[2, end] / wf_right[2, 1]
-        @assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
-        isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
-        wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
-        wf = hcat(wf[:, 1:(end - 1)], wf_right_rescaled)
-    end
+    # left partition
+    prob = ODEProblem(dirac!, [0, 1], (0, r_mid))
+    sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r)
+    wf_left = hcat(sol.u...)
+
+    # right partition
+    prob = ODEProblem(dirac!, [0, 1], (s.r_mesh.r_max, r_mid))
+    sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
+    wf_right = reverse(hcat(sol.u...); dims=2)    
+
+    # join two segments
+    rescale_factor_g = wf_left[1, end] / wf_right[1, 1]
+    rescale_factor_f = wf_left[2, end] / wf_right[2, 1]
+    @assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
+    isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
+    wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
+    wf = hcat(wf_left[:, 1:(end - 1)], wf_right_rescaled)
 
     if normalize
-        norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
+        norm = sum(wf .* wf .* transpose(s.r_mesh.w)) # integration by Reimann sum
         wf = wf ./ sqrt(norm)
     end
 
@@ -80,10 +81,10 @@ end
 "Returns a function that solves the Dirac equation in two partitions and returns 
     the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
     where is r is in fm."
-function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Tsit5())
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_mesh.r_max/2, algo=Tsit5())
     (f1, f2) = optimized_dirac_potentials(p, s)
     prob_left  = ODEProblem(dirac!, [0.0, 1.0], (0, r))
-    prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_max, r))
+    prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_mesh.r_max, r))
     function func(E)
         u_left  = solve(prob_left,  algo, p=(κ, E, f1, f2), saveat=[r])
         u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r])
@@ -141,22 +142,22 @@ j_κ(κ::Int) = abs(κ) - 1/2
 "Orbital angular momentum l for a given κ value"
 l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 
-"Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
+"Calculate scalar and vector densities of a nucleon species evaluated at mesh points and returns them as vectors (ρ_s, ρ_v) where
     the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
     the other parameters are defined above"
 function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}}
     spectrum = p ? s.p_spectrum : s.n_spectrum
     (κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
     
-    ρr2 = zeros(2, s.divs + 1) # ρ×r² values
+    ρr2 = zeros(2, length(s.r_mesh.r)) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
-        wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
+        wf = solveNucleonWf(κ, p, E, s; normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end
 
-    r2s = rs(s).^2
+    r2s = s.r_mesh.r .^ 2
     ρ = ρr2 ./ transpose(r2s)
     ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
     
diff --git a/system.jl b/system.jl
index d52f6c2..39d916b 100644
--- a/system.jl
+++ b/system.jl
@@ -10,13 +10,26 @@ end
 "Initializes an unfilled spectrum"
 unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
 
+"Defines a mesh"
+struct mesh
+    r_max::Float64
+    r::Vector{Float64}
+    w::Vector{Float64}
+end
+
+"Create a uniform mesh"
+function uniform_mesh(r_max::Float64, divs::Int)::mesh
+    r = range(0, r_max, length=divs+1) |> collect
+    w = fill(r_max / divs, divs+1)
+    return mesh(r_max, r, w)
+end
+
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system
     Z::Int
     N::Int
 
-    r_max::Float64
-    divs::Int
+    r_mesh::mesh
 
     p_spectrum::spectrum
     n_spectrum::spectrum
@@ -32,7 +45,7 @@ mutable struct system
     ρ_vn::Vector{Float64}
 
     "Initialize an unsolved system"
-    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
+    system(Z, N, r_max, divs) = new(Z, N, uniform_mesh(r_max, divs), unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
 
     "Dummy struct to define the mesh"
     system(r_max, divs) = system(0, 0, r_max, divs)
@@ -41,14 +54,8 @@ end
 "Get mass number of nucleus"
 A(s::system)::Int = s.Z + s.N
 
-"Get r values in the mesh"
-rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)
-
-"Get Δr value for the mesh"
-Δr(s::system)::Float64 = s.r_max / s.divs
-
 "Get the number of protons or neutrons in the system"
 Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N
 
 "Create an empty array for the size of the mesh"
-zero_array(s::system) = zeros(1 + s.divs)
+zero_array(s::system) = zeros(length(s.r_mesh.r))
diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl
index 1e11111..5dc335a 100644
--- a/test/Pb208_nucleon_dens.jl
+++ b/test/Pb208_nucleon_dens.jl
@@ -23,13 +23,13 @@ s.A0 = As
 
 solveNucleonDensity!(true, s)
 
-p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+p_sp = plot(s.r_mesh.r, s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_vp = plot(s.r_mesh.r, s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
 solveNucleonDensity!(false, s)
 
-p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_sn = plot(s.r_mesh.r, s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_vn = plot(s.r_mesh.r, s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
 
 # benchmark data generated from Hartree.f
 # format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)
diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl
index 70aa5ae..fa58906 100644
--- a/test/Pb208_nucleon_wf.jl
+++ b/test/Pb208_nucleon_wf.jl
@@ -28,6 +28,6 @@ wf = solveNucleonWf(κ, p, groundE, s)
 gs = wf[1, :]
 fs = wf[2, :]
 
-plot(rs(s), gs, label="g(r)")
-plot!(rs(s), fs, label="f(r)")
+plot(s.r_mesh.r, gs, label="g(r)")
+plot!(s.r_mesh.r, fs, label="f(r)")
 xlabel!("r (fm)")