(*The product of two 1s orbitals. Note that the 1sa1sb density represents
 a kind of rope of nearly equal thickness along the bond, with an exponential
 decay perpendicularly to the bond. Outside the bond region the density of this rope 
decays very fast.According to the analysis p.English387 the chemical
bond stems from the interaction of the nuclear charge with this rope. Thus,
the popular representation of the chemical bond as a stick connecting the atomic balls
is not as far from reality as one might think*)
R=4.
pi=3.141592532 
a=(Sqrt[(1/pi)]) (E^(-(Sqrt[(x^2)+(y^2)])))
b=(Sqrt[(1/pi)]) (E^(-(Sqrt[(((x-R))^2)+(y^2)])))
Plot3D[a*b,{x,-2.,6.},{y,-2.,2.},PlotPoints->24]
Plot3D[a^2+b^2,{x,-2.,6.},{y,-2.,2.},PlotPoints->24]
Plot3D[a^2+2*a*b+b^2,{x,-2.,6.},{y,-2.,2.},PlotPoints->24]
Plot3D[a^2-2*a*b+b^2,{x,-2.,6.},{y,-2.,2.},PlotPoints->24]

R=2.
pi=3.141592532 
a=(Sqrt[(1/pi)]) (E^(-(Sqrt[(x^2)+(y^2)])))
b=(Sqrt[(1/pi)]) (E^(-(Sqrt[(((x-R))^2)+(y^2)])))
Plot3D[a*b,{x,-2.,6.},{y,-2.,2.},PlotPoints->24]
R=1.
pi=3.141592532 
a=(Sqrt[(1/pi)]) (E^(-(Sqrt[(x^2)+(y^2)])))
b=(Sqrt[(1/pi)]) (E^(-(Sqrt[(((x-R))^2)+(y^2)])))
Plot3D[a*b,{x,-2.,6.},{y,-2.,2.},PlotPoints->24,PlotRange->{0.,.12}]
R=0.1
pi=3.141592532 
a=(Sqrt[(1/pi)]) (E^(-(Sqrt[(x^2)+(y^2)])))
b=(Sqrt[(1/pi)]) (E^(-(Sqrt[(((x-R))^2)+(y^2)])))
Plot3D[a*b,{x,-2.,6.},{y,-2.,2.},PlotPoints->24,PlotRange->{0.,.30}]