diff --git a/B2R_comparison.jl b/B2R_comparison.jl
new file mode 100644
index 0000000..3b12791
--- /dev/null
+++ b/B2R_comparison.jl
@@ -0,0 +1,51 @@
+using Plots, LinearAlgebra
+include("p_space.jl")
+
+berggren_mesh = get_mesh((0, 0.4 - 0.15im, 0.8, 6), 128)
+csm_mesh = get_mesh((0, 8 - 3im), 512)
+
+for mesh in (berggren_mesh, csm_mesh)
+    p, w = mesh
+    mesh_E = p.*p ./ (2*0.5)
+
+    # ResonanceEC: Eq. (20)
+    V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
+
+    training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
+    training_E = Vector{ComplexF64}(undef, length(training_points))
+    EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
+
+    for (j, c) in enumerate(training_points)
+        evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
+        i = identify_pole_i(p, evals)
+        training_E[j] = evals[i]
+        EC_basis[:, j] = evecs[:, i]
+    end
+
+    EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
+    EC_basis_w = EC_basis .* w
+    N_EC = transpose(EC_basis_w) * EC_basis
+
+    extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+
+    exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
+    extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
+
+    for (j, c) in enumerate(extrapolate_points)
+        exact_E[j] = quick_pole_E(V_system(c))
+
+        H = get_H_matrix(V_system(c), p, w)
+        H_EC = transpose(EC_basis_w) * H * EC_basis
+        evals = eigvals(H_EC, N_EC)
+        i = argmin(abs.(evals .- exact_E[j]))
+        extrapolate_E[j] = evals[i]
+    end
+
+    scatter(real.(training_E), imag.(training_E), label="training")
+    scatter!(real.(exact_E), imag.(exact_E), label="exact")
+    scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
+    plot!(real.(mesh_E), imag.(mesh_E), label="contour")
+    xlims!(-0.3,0.3)
+    ylims!(-0.120,0.020)
+    savefig("temp/" * string(rand(UInt16)) * ".pdf")
+end
\ No newline at end of file
diff --git a/B2R_test.jl b/B2R_test.jl
deleted file mode 100644
index 111ac1f..0000000
--- a/B2R_test.jl
+++ /dev/null
@@ -1,46 +0,0 @@
-using Plots, LinearAlgebra
-include("p_space.jl")
-
-vertices = (0, 0.4 - 0.15im, 0.8, 6)
-subdivisions = 128
-p, w = get_mesh(vertices, subdivisions)
-mesh_E = p.*p ./ (2*0.5)
-
-# ResonanceEC: Eq. (20)
-V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
-
-training_points = range(1.35, 0.9, 5)
-training_E = Vector{ComplexF64}(undef, length(training_points))
-EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
-
-for (j, c) in enumerate(training_points)
-    evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
-    i = identify_pole_i(p, evals)
-    training_E[j] = evals[i]
-    EC_basis[:, j] = evecs[:, i]
-end
-
-EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
-EC_basis_w = EC_basis .* w
-N_EC = transpose(EC_basis_w) * EC_basis
-
-extrapolate_points = range(0.75, 0.40, 8)
-
-exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
-extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
-
-for (j, c) in enumerate(extrapolate_points)
-    exact_E[j] = quick_pole_E(V_system(c))
-
-    H = get_H_matrix(V_system(c), p, w)
-    H_EC = transpose(EC_basis_w) * H * EC_basis
-    evals = eigvals(H_EC, N_EC)
-    i = argmin(abs.(evals .- exact_E[j]))
-    extrapolate_E[j] = evals[i]
-end
-
-scatter(real.(training_E), imag.(training_E), label="training")
-scatter!(real.(exact_E), imag.(exact_E), label="exact")
-scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
-plot!(real.(mesh_E), imag.(mesh_E), label="contour")
-xlims!(-0.7,0.7)
\ No newline at end of file