In[1]:= Psi1[l1_,m1_,l2_,m2_]:=Sqrt[(2*l1+1)/(4*Pi)]*LegendreP[l1,m1,Cos[theta1]]*Sqrt[(2*l2+1)/(4*Pi)]*LegendreP[l2,m2,Cos[theta2]]-Sqrt[(2*l1+1)/(4*Pi)]*LegendreP[l1,m1,Cos[theta2]]*Sqrt[(2*l2+1)/(4*Pi)]*LegendreP[l2,m2,Cos[theta1]]

In[2]:= Den1[l1_,m1_,l2_,m2_]:=1/2*Psi1[l1,m1,l2,m2]^2

In[3]:= Dens1[l1_,m1_,l2_,m2_]:=Plot3D[Integrate[Den1[l1,m1,l2,m2],{phi1,0,2*Pi},{phi2,0,2*Pi}]*Sin[theta1]*Sin[theta2],{theta1,0,Pi},{theta2,0,Pi}]

In[4]:= Do[Print[{l1,l2,Dens1[l1,m1,l2,m2]}],{l1,0,5},{l2,0,5},{m1,0,0},{m2,0,0}]