Vtot[x_,y_,z_]=1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]

ContourPlot3D[Vtot[x,y,z],{x,-3,3},{y,-3,3},{z,-3,3}]

ContourPlot3D[1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]==0,{x,-2,2},{y,-2,2},{z,-2,2}]
ContourPlot3D[1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]==.05,{x,-2,2},{y,-2,2},{z,-2,2}]
ContourPlot3D[1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]==.1,{x,-2,2},{y,-2,2},{z,-2,2}]
ContourPlot3D[1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]==.2,{x,-2,2},{y,-2,2},{z,-2,2}]
ContourPlot3D[1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]==.5,{x,-2,2},{y,-2,2},{z,-2,2}]
ContourPlot3D[1/Sqrt[x^2+y^2+(z-1)^2] - 1/Sqrt[x^2+y^2+(z+1)^2]==1,{x,-2,2},{y,-2,2},{z,-2,2}]

Ez[x_,y_,z_]=(z-1)/(x^2+y^2+(z-1)^2)^(3/2) - (z+1)/(x^2+y^2+(z+1)^2)^(3/2)
Ezonaxis[z_]= Ez[x,y,z]/. {x->0,y->0}

((-1+z)/((((((-1+z))^2)))^(3/2)))-((1+z)/((((((1+z))^2)))^(3/2)))
Plot[Ezonaxis[z],{z,0,4},PlotRange-> {-100,100}]
Plot[Ezonaxis[z],{z,0,4},PlotRange-> {-10,10}]