diff --git a/XZ_technique.jl b/XZ_technique.jl
index b037a49..284ca69 100644
--- a/XZ_technique.jl
+++ b/XZ_technique.jl
@@ -8,7 +8,7 @@ V1 = -8.6
 R1 = 1.75
 V2 = 3.4
 R2 = 3.2
-n_max = 20
+n_max = 30
 
 ns = collect(0:n_max)
 ls = fill(l, n_max + 1)
diff --git a/ho_basis.jl b/ho_basis.jl
index 0074a0b..5948687 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -2,10 +2,11 @@ using NuclearToolkit
 using SpecialFunctions
 
 # Gaussian potentials in HO space
+inv_factorial(n) = Iterators.prod(inv.(1:n))
 sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
 sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
 N_nl(n, l) = (-1)^n * sqrt_factorial(n) * sqrt(1/sqrt(pi) * 2^(n+l+2)) / sqrt_double_factorial(2*n+2*l+1)
-prefactor(n, l, k) = (-1)^k * binomial(n + l + 1/2, n - k) / factorial(k)
+prefactor(n, l, k) = (-1)^k * binomial(n + l + 1/2, n - k) * inv_factorial(k)
 Talmi(l, R, k1, k2) = (1/2) / (1 + 1/R^2)^(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
 V_Gaussian(R, l, n1, n2) = N_nl(n1, l) * N_nl(n2, l) * sum([prefactor(n1, l, k1) * prefactor(n2, l, k2) * Talmi(l, R, k1, k2) for (k1, k2) in Iterators.product(0:n1, 0:n2)])