Clear["Global`*"]

H0 - Matrix Elements

Clear[l];
f[r_]:=r^l Exp[-bj r^2];
g[r_]:=D[f[r],r];
h[r_]:=D[g[r],r];
Table[Simplify[Assuming[Re[bi+bj]>0 && Re[l]≥0,Integrate[r^(l+2) Exp[-bi r^2] (-h[r]/2-2/r g[r]/2+l(l+1)/r^2 /2f[r]-Z/r f[r]),{r,0,∞}]]],{l,0,2}]

S - Matrix Elements

Table[Assuming[Re[bi+bj]>0 && Re[l]≥0,Integrate[r^(2l+2)Exp[-bi r^2]Exp[-bj r^2],{r,0,∞}]],{l,0,2}]

J - Matrix Elements

Clear[l1,l2,m1,m2,p,q,JF,JG,JH];
JF=Assuming[Re[q]>0 && Re[p]>0&& Re[l]≥0 && Re[l1]≥0 && Re[l2]≥0, Integrate[r1^(2l1+2)Exp[-p r1^2](Integrate[r2^(2l2+2)Exp[-q r2^2] ((r2^l)/(r1^(l+1))),{r2,0,r1}]+Integrate[r2^(2l2+2)Exp[-q r2^2] ((r1^l)/(r2^(l+1))),{r2,r1,∞}]),{r1,0,∞}]];

JG:=
(1/(2l1+1))*
Sum[Integrate[Integrate[Conjugate[SphericalHarmonicY[l1,m1,t,p]]SphericalHarmonicY[l,0,t,p]SphericalHarmonicY[l1,m1,t,p] Sin[t],{t,0,Pi}],{p,0,2Pi}],{m1,-l1,l1}];
JH:=(1/(2l2+1))*Sum[Integrate[Integrate[Conjugate[SphericalHarmonicY[l2,m2,t,p]]SphericalHarmonicY[l,0,t,p]SphericalHarmonicY[l2,m2,t,p] Sin[t],{t,0,Pi}],{p,0,2Pi}],{m2,-l2,l2}];

Table[Table[Simplify[Sum[((4Pi)/(2l+1))JF*JG*JH,{l,0,l1+l2}],Assumptions->Re[p]>0&&Re[q]>0 &&Im[q]==0&&Im[p]==0],{l1,0,1}],{l2,0,1}]

K - Matrix Elements

Clear[l1,l2,m1,m2,p,q,KF,KG,KH];
KF=Assuming[Re[q]>0 && Re[p]>0&& Re[l]>=0 && Re[l1]≥0 && Re[l2]≥0 , Integrate[r1^(l1+l2+2)Exp[-p r1^2](Integrate[r2^(l1+l2+2)Exp[-q r2^2] ((r2^l)/(r1^(l+1))),{r2,0,r1}]+Integrate[r2^(l1+l2+2)Exp[-q r2^2] ((r1^l)/(r2^(l+1))),{r2,r1,∞}]),{r1,0,∞}]];

KG:=1/Sqrt[2l2+1]*Sum[Integrate[Integrate[Conjugate[SphericalHarmonicY[l1,m1,t,p]]SphericalHarmonicY[l,m2-m1,t,p]SphericalHarmonicY[l2,m2,t,p] Sin[t],{t,0,Pi}],{p,0,2Pi}],{m2,-l2,l2}];
KH:=1/Sqrt[2l2+1]*Sum[Integrate[Integrate[Conjugate[SphericalHarmonicY[l2,m2,t,p]]Conjugate[SphericalHarmonicY[l,m2-m1,t,p]]SphericalHarmonicY[l1,m1,t,p] Sin[t],{t,0,Pi}],{p,0,2Pi}],{m2,-l2,l2}];

Clear[m1,m2,l1,l2]
m1=0;
Table[Table[Simplify[Sum[((4Pi)/(2l+1))KF *Conjugate[KG]KG,{l,0,l1+l2}],Assumptions->Re[p]>0&&Re[q]>0 &&Im[q]==0&&Im[p]==0],{l1,0,1}],{l2,0,1}]

l1=0;l2=1;
Simplify[Sum[((4Pi)/(2l+1))KH KG,{l,0,l1+l2}],Assumptions->Re[p]>0&&Re[q]>0 &&Im[q]==0&&Im[p]==0]
l1=1; l2=0;
Simplify[Sum[((4Pi)/(2l+1))KH KG,{l,0,l1+l2}],Assumptions->Re[p]>0&&Re[q]>0 &&Im[q]==0&&Im[p]==0]