Help with Mathematica for optimizing a Gaussian basis set for a hydrogen atom in a magnetic field

Question

This is a follow-up to my previous question: Optimal Gaussian basis set for hydrogen atom in magnetic field

Brief description of the problem

I would like to find the ground and excited states of the system:

$$H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{B^2}{8}\rho^2-\frac{B}{2} (m + 2 m_s),\tag{1}$$ in which $r=(\rho,z,\phi), B=5, m=0, m_s=-1/2$.

Method

To find the eigenvalues I would like to use Gaussian functions: $$\psi_{ij}=e^{-b_{j} z^2}e^{-a_{i} \rho^2},\tag{2}$$

in which,

\begin{align} a_{i}&=a_1 qa^{i-1}, i=1, 2 ,3,..., i_{\textrm{max}},\tag{3}\\b_{j}&=b_1 qb^{j-1}, j=1, 2 ,3,..., j_{\textrm{max}},\tag{4} \end{align}

are geometrical progressions.

Code

I used Mathematica (in the code I rename $\alpha_i \equiv a_i$, $\beta_j \equiv b_j$, $\rho\equiv r$):

ClearAll["Global`*"]

imax = 4; jmax = 5;

B = 5; lz = m = 0; ms = -1/2;

VP1[r_, z_] := -1/Sqrt[r^2 + z^2];

Psi[r_, z_, i_, j_] := Exp[-b[j]*z^2]*Exp[-a[i]*r^2];

a[1] = a1;
b[1] = b1;
Do[a[i] = a[i - 1] qa, {i, 2, imax}];
Do[b[j] = b[j - 1] qb, {j, 2, jmax}];

(*kinetic energy*)
Kk = FullSimplify[
   Psi[r, z, i2, j2] *
    Laplacian[Psi[r, z, i1, j1] , {r, \[Theta], z}, "Cylindrical"]];
Kk1 = -1/2 2 Pi*
   Integrate[Kk r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
      a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
Kx = Table[ 
   Kk1, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];
KK[a1_, b1_, qa_, qb_] = Flatten[Kx, {{1, 3}, {2, 4}}];

(*potential energy + A kinetic energy*)
Px1 = 2 Pi*
   Integrate[
    Psi[r, z, i2, j2] *(VP1[r, z] + 1/8 B^2 r^2 + B/2 (m + 2 ms))*
     Psi[r, z, i1, j1]*r, {r, 0, \[Infinity]}, {z, -Infinity, 
     Infinity}, 
    Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
      a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
Px = Table[
   Px1, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];
PP = Flatten[Px, {{1, 3}, {2, 4}}];

EE = KK + PP;

(*normalization*)
int = 2 Pi*
   Integrate[
    Psi[r, z, i2, j2] *Psi[r, z, i1, j1]  r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}, 
    Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
      a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
norm = Table[ 
   int, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}];
BB = Flatten[norm, {{1, 3}, {2, 4}}];

IBB = Inverse[BB](**)

(*the sum of eigenvalues to minimize*)
Tr[EE.IBB]

Problem

The code cannot handle more than 20 basis functions, but the paper from the 1970s was able to do 100+ basis functions. I will be glad to see any comments. Please use any programming language convenient for you, I will try to figure it out.