Z = 1
r = (Sqrt[(x^2) + (y^2)])
cz = (E^(-((Z r)/2)))
Ns = (Sqrt[Z^(3/(32 π))])
Np = (Sqrt[Z^(5/(32 π))])
px = Np x cz
s = Ns (2 - Z r) cz
u = 0
h =  s + u px
ContourPlot[h, {x, -18, 18}, {y, -18, 18}, 
      PlotPoints -> 60]
u = 1410
h =  s + u px
ContourPlot[h, {x, -18, 18}, {y, -18, 18}, 
      PlotPoints -> 60]
u = Sqrt[3]
hsp3 =  s + u px
ContourPlot[hsp3, {x, -18, 18}, {y, -18, 18}, 
      PlotPoints -> 60]
u = Sqrt[2]
hsp2 = s + u px
ContourPlot[hsp2, {x, -18, 18}, {y, -18, 18}, 
      PlotPoints -> 60]
u = Sqrt[1]
hsp = s + u px
ContourPlot[hsp, {x, -18, 18}, {y, -18, 18}, 
      PlotPoints -> 60]