R=2 
Clear[x, y]
r = (Sqrt[(x^2) + (y^2)])
A = (Sqrt[(1/π)]) (E^(-r))
B = (Sqrt[(1/π)]) (E^(-(Sqrt[(((x - R))^2) + (y^2)])))
S = (E^(-R)) (((R^2)/3) + R + 1)
Fi1 = ((A + B)/(Sqrt[2 (1 + S)]))
Fi2 = ((A - B)/(Sqrt[2 (1 - S)]))
Psi1 = (((((Sqrt[(1/π)]))^2) (Fi2 - Fi1))/(Sqrt[6]))
Psikw1 = (Psi1^2)
Psi2 = (((((Sqrt[(1/π)]))^2) (Fi1 + Fi2))/(Sqrt[6]))
Psikw2 = (Psi2^2)
Plot3D[Psikw1, {x, -1, 3}, {y, -2, 2}, 
      PlotPoints -> 50]
Plot3D[Psikw2, {x, -1, 3}, {y, -2, 2}, 
      PlotPoints -> 50]