diff --git a/3body_resonance.jl b/3body_resonance.jl
index 851d677..248eae1 100644
--- a/3body_resonance.jl
+++ b/3body_resonance.jl
@@ -1,4 +1,4 @@
-using Arpack, SparseArrays
+using Arpack, SparseArrays, LRUCache
 include("ho_basis.jl")
 include("p_space.jl")
 
@@ -38,8 +38,8 @@ atol = 10^-6
 maxevals = 10^5
 V1_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μ1ω1, atol=atol, maxevals=maxevals)
 V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
-V1_cache = fill(Complex(NaN), 1+l_max, 1+n_max, 1+n_max)
-V_relative_cache = fill(Complex(NaN), 1+l_max, 1+n_max, 1+n_max)
+V1_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
+V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
 
 @time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V1_cache)
 @time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_relative_cache) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2, dtype=ComplexF64, cache=V_relative_cache)
diff --git a/ho_basis.jl b/ho_basis.jl
index c22431c..84dda6f 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,6 +1,7 @@
 using SparseArrays
 using SpecialFunctions
 using QuadGK
+using LRUCache
 include("helper.jl")
 
 # Gaussian potentials in HO space
@@ -71,19 +72,15 @@ function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), μω_gen=1.0, μ
     return (μω_gen / μ) .* mat
 end
 
-function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64, cache=fill(convert(dtype, NaN), 1+maximum(ls), 1+maximum(ns), 1+maximum(ns)))
+function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64, cache=LRU{Tuple{UInt8, UInt8, UInt8}, dtype}(maxsize=(1+maximum(ns))^2))
     mat = zeros(dtype, length(ns), length(ns))
     Threads.@threads for idx in CartesianIndices(mat) 
         if !mask[idx]; continue; end
         (i, j) = Tuple(idx)
         if ls[i] == ls[j]
-            l = ls[i]
-            n1, n2 = minmax(ns[i], ns[j]) # assuming transpose symmetry
-            if isnan(cache[1+l, 1+n1, 1+n2])
-                cache[1+l, 1+n1, 1+n2] = V_l(l, n1, n2) # hopefully no race condition
-                @assert !isnan(cache[1+l, 1+n1, 1+n2]) "V matrix element returned NaN"
-            end
-            mat[idx] = cache[1+l, 1+n1, 1+n2]
+            l = UInt8(ls[i])
+            n1, n2 = UInt8.(minmax(ns[i], ns[j])) # assuming transpose symmetry
+            mat[idx] = (get!(cache, (l, n1, n2)) do; V_l(l, n1, n2); end)
         end
     end
     return sparse(mat)
diff --git a/ho_basis_3body.jl b/ho_basis_3body.jl
index 47882a4..d9f57cd 100644
--- a/ho_basis_3body.jl
+++ b/ho_basis_3body.jl
@@ -1,4 +1,4 @@
-using Arpack, SparseArrays
+using Arpack, SparseArrays, LRUCache
 include("ho_basis.jl")
 include("p_space.jl")
 
@@ -35,8 +35,8 @@ println("Constructing KE matrices")
 println("Constructing PE matrices")
 V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μ1ω1)
 V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_global)
-V1_cache = fill(NaN, 1+l_max, 1+n_max, 1+n_max)
-V_relative_cache = fill(NaN, 1+l_max, 1+n_max, 1+n_max)
+V1_cache = LRU{Tuple{UInt8, UInt8, UInt8}, Float64}(maxsize=(1+l_max)*(1+n_max)^2)
+V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, Float64}(maxsize=(1+l_max)*(1+n_max)^2)
 @time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s; mask=mask1, cache=V1_cache)
 @time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, cache=V_relative_cache) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2, cache=V_relative_cache)
 @time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)