Two electron integrals between 1s STO
@TAR86 was correct, I mistook the volume integral as simple 1d integral. Also integrating $1/|r_2 - r_1|$ would result in a singularity at $r_1 == r_2$ which would yield really wrong results.
As per @user1271772 suggestions, materials modelling stack exchange was quite useful (see question here). There as per Susi Lehota's suggestion I was able to get correct result by integrating Legendre Expansion.
My octave code for two electron integral (11|11) follows:
zeta1=1.45;
chi1111 = @(r1,r2) 16*zeta1^3*zeta1^3*exp(-2*zeta1.*(r1)).*exp(-2*zeta1.*(r2)).*r1.*r1.*r2.*r2./max(r1,r2);
dblquad(chi1111,0.000001,10,0.0000001,10,1E-3)
ans =  0.90623

Similarly other integrals can be encoded, and I get results exactly as Levine example.
Thank you everyone