diff --git a/ho_basis_3body.jl b/ho_basis_3body.jl
new file mode 100644
index 0000000..5cac534
--- /dev/null
+++ b/ho_basis_3body.jl
@@ -0,0 +1,41 @@
+using Arpack, SparseArrays
+include("ho_basis.jl")
+include("p_space.jl")
+
+E_max = 10
+ω = 1.0
+Λ = 0
+m = 1.0
+
+Va = -5
+Ra = sqrt(3)
+Vb = 2
+Rb = sqrt(10)
+
+μ1 = m * 1/2
+μ2 = m * 2/3
+c1 = 1/sqrt(2)
+c2 = sqrt(2)
+
+Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
+println("Basis size = ", length(Es))
+
+println("Constructing KE matrices")
+T1 = sp_T_matrix(n1s, l1s; ω=ω, μ=c1^2 * μ1)
+T2 = sp_T_matrix(n2s, l2s; ω=ω, μ=c2^2 * μ2)
+
+println("Constructing PE matrices")
+V_l(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω) + Vb * V_Gaussian(Rb, l, n1, n2; ω=ω)
+V1 = sp_V_matrix(V_l, n1s, l1s)
+V_relative = sp_V_matrix(V_l, n1s, l1s) + sp_V_matrix(V_l, n2s, l2s)
+
+println("Performing Moshinsky transform")
+@time "Calculating brackets" Moshinsky_mat = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
+@time "Applying transform" V2 = Moshinsky_mat * V_relative
+
+H = T1 + T2 + V1 + V2
+
+println("Diagonalizing H matrix")
+@time "Diagonalizing" evals, _ = eigs(H, nev=5, which=:SM)
+
+display(evals)
\ No newline at end of file