diff --git a/ho_basis.jl b/ho_basis.jl
index 82d8238..3243755 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -54,8 +54,6 @@ function get_2p_basis(E_max)
     return (Es, n1s, l1s, n2s, l2s)
 end
 
-get_V_matrix(V_l, ls, ns) = throw("unimplemented")
-
 function sp_T_matrix(ns, ls; ω=1.0, μ=1.0)
     mat = spzeros(length(ns), length(ns))
     for idx in CartesianIndices(mat)
@@ -74,6 +72,17 @@ function sp_T_matrix(ns, ls; ω=1.0, μ=1.0)
     return (ω / μ) .* mat
 end
 
+function sp_V_matrix(V_l, ns, ls; dtype=Float64)
+    mat = spzeros(dtype, length(ns), length(ns))
+    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
+end
+
 get_H_matrix(V_l, ns, ls) = get_T_matrix(ns, ls) + get_V_matrix(V_l, ns, ls)
 
 function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
diff --git a/ho_basis_benchmark.jl b/ho_basis_benchmark.jl
index ab2c1e8..58c3309 100644
--- a/ho_basis_benchmark.jl
+++ b/ho_basis_benchmark.jl
@@ -15,7 +15,9 @@ ns = collect(0:n_max)
 ls = fill(l, n_max + 1)
 
 T = sp_T_matrix(ns, ls; ω=ω, μ=μ)
-V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω)
+
+V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; ω=ω) + V2 * V_Gaussian(R2, l, n1, n2; ω=ω)
+V = sp_V_matrix(V_l, ns, ls; dtype=ComplexF64)
 
 cs = range(1.25, 0.25, 10)
 
@@ -24,7 +26,7 @@ bench_E = similar(cs, ComplexF64)
 
 for (j, c) in enumerate(cs)
     H = T + c .* V
-    evals = eigvals(H)
+    evals = eigvals(collect(H))
     bench_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.4, meshpoints=512)
     i = argmin(abs.(evals .- bench_E[j]))
     E[j] = evals[i]
diff --git a/ho_basis_test.jl b/ho_basis_test.jl
index 439f4c9..6d26fb0 100644
--- a/ho_basis_test.jl
+++ b/ho_basis_test.jl
@@ -10,7 +10,10 @@ ns = collect(0:n_max)
 ls = fill(l, n_max + 1)
 
 T = sp_T_matrix(ns, ls)
-V = V0 .* V_Gaussian.(R, l, ns, transpose(ns))
+
+V_l(l, n1, n2) = V0 * V_Gaussian(R, l, n1, n2)
+V = sp_V_matrix(V_l, ns, ls)
+
 H = T + V
 
-eigvals(H)
\ No newline at end of file
+eigvals(collect(H))
\ No newline at end of file