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)])))
S = (1 + R + ((R^2)/3)) (E^(-R))
Heitler = (((Sqrt[(1/pi)]) (b + (E^(-R)) a))/(Sqrt[2 (1 + (S^2))]))
Plot3D[Heitler, {x, -2. , 6. }, {y, -2. , 2. }, 
      PlotPoints -> 25]
HFock = (((Sqrt[(1/pi)]) (a + b + (E^(-R)) a + (E^(-R)) b))/((Sqrt[1 + (S^2)]) 2))
Plot3D[HFock, {x, -2. , 6. }, {y, -2. , 2. }, 
      PlotPoints -> 25]