diff --git a/berggren_2body.jl b/berggren_2body.jl
new file mode 100644
index 0000000..95017be
--- /dev/null
+++ b/berggren_2body.jl
@@ -0,0 +1,25 @@
+using SparseArrays, Arpack
+include("helper.jl")
+include("p_space.jl")
+
+E_target = -0.3919
+
+μ = 0.5
+Va = -2
+Ra = 2
+V_of_r(r) = Va * exp(-r^2 / Ra^2)
+
+vertices = [0, 0.5 - 0.3im, 1, 4]
+subdivisions = [16, 16, 16]
+ks, ws = get_mesh(vertices, subdivisions)
+ls = collect(0:4)
+
+V_l(l, k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)
+
+H_l = [get_H_matrix((k, kp) -> V_l(l, k, kp), ks, ws, μ) for l in ls]
+
+H1 = blockdiag(sparse.(H_l)...)
+
+H = H1
+@time "Eigenvalues" evals, _ = eigs(H, nev=10, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+E = nearest(evals, E_target)