Overlap matrix (H2 molecule, STO-3G)
You have two errors in your code, one minor and one major.
The minor error is in your multiplication function. The value you return should have $\pi^{3/2}$, but the way it is arranged now the power is only being applied to the denominator of that expression.
The major error is how you are normalizing your overlap integrals. You can't simply divide $\langle1|2\rangle$ by $\langle1|1\rangle$ in order to normalize it. You have to explicitly normalize each of the Gaussian basis functions. To do this you can include a factor of ((4*self.alpha*other.alpha)/np.pi**2)**(3/4) in the value you are returning from your multiplication (in terms of readability of your code, it may be wiser to define variables within the function like proportionality_constant, gauss_integral, and normalization and then return the product of these). Or you could simplify and just multiply the exponential by ((4*self.alpha*other.alpha)/(self.alpha+other.alpha)**2)**(3/4) This also means you only need the second part of your overlap function, since the multiplication will take care of the normalization.
The reason your code works for STO-1G is that dividing by $\langle1|1\rangle$ accidentally gives the correct normalization since there is only one normalization constant associated with the one Gaussian function of STO-1G.