diff --git a/ho_basis.jl b/ho_basis.jl
index 0d1bcf0..9b087e4 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,21 +1,8 @@
 using SparseArrays
-using SpecialFunctions
 using QuadGK
 using LRUCache
 include("helper.jl")
-
-# Gaussian potentials in HO space
-inv_factorial(n) = Iterators.prod(inv.(1:n))
-sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
-N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
-Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
-
-# numerical evaluation of V matrix elements
-sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
-sqrt_sqrt_pi = sqrt(sqrt(pi))
-laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k)  for k in 0:n])
-ho_basis(l, n, x) = (-1)^n / sqrt_sqrt_pi * 2^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1) * x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
+include("math.jl")
 
 function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
     integrand(r) = sqrt(μω_gen) * ho_basis(l, n1, sqrt(μω_gen) * r) * ho_basis(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
@@ -23,8 +10,6 @@ function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
     return integral
 end
 
-##############################################################
-
 function get_sp_basis(E_max)
     Es = Int[]
     ns = Int[]
diff --git a/math.jl b/math.jl
new file mode 100644
index 0000000..331af70
--- /dev/null
+++ b/math.jl
@@ -0,0 +1,14 @@
+using SpecialFunctions
+
+# Gaussian potentials in HO space
+inv_factorial(n) = Iterators.prod(inv.(1:n))
+sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
+N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
+Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
+V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+
+# for numerical evaluation of V matrix elements
+sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
+sqrt_sqrt_pi = sqrt(sqrt(pi))
+laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k)  for k in 0:n])
+ho_basis(l, n, x) = (-1)^n / sqrt_sqrt_pi * 2^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1) * x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)