diff --git a/ho_basis.jl b/ho_basis.jl
index a5f44e5..8547e5a 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -71,14 +71,14 @@ function sp_T_matrix(ns, ls; ω=1.0, μ=1.0)
 end
 
 function sp_V_matrix(V_l, ns, ls; dtype=Float64)
-    mat = spzeros(dtype, length(ns), length(ns))
-    for idx in CartesianIndices(mat)
+    mat = zeros(dtype, length(ns), length(ns))
+    Threads.@threads for idx in CartesianIndices(mat)
         (i, j) = Tuple(idx)
         if ls[i] == ls[j]
             mat[idx] = V_l(ls[i], ns[i], ns[j])
         end
     end
-    return mat
+    return sparse(mat)
 end
 
 function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)