<<"VectorAnalysis`"

" We define the function ";
phi[r_,ri_]:=(q/(Sqrt[DotProduct[r-ri,r-ri]]))

" this can be used to define any potential function due   to a point charge at ri ";

potential1=phi[{Xx,Yy,Zz},{0,0,0}]/.q->1

ContourPlot[potential1/.Zz->0.0,{Xx,-1,1},{Yy,-1,1},ContourShading->False,PlotPoints->40]

efield1=-Grad[potential1]

<<"VectorFieldPlots`"

Needs["VectorFieldPlots`"];

VectorFieldPlots`VectorFieldPlot[{efield1[[1]],efield1[[2]]}/.Zz->0,{Xx,-1,1},{Yy,-1,1},ScaleFactor->0.1`]

Some fun things

Linear Superposition; use Mathematica to 
1) Find an expression for the potential due to  a charge q=+1, at position r={0.5,0,0}.

2) Find an expression for the potential due to a charge q=-1, at position r={-0.5,0,0}.

3) Construct the electric field for each part 1) 2)

4) Add the results of 1) and 2) to find the potential  for a dipole configuration.

5) Do a contour plot of this potential, on the z=0 plane

6) Find the electric field for potential 5), and plot it on the z=0 plane.

7) Use the discussion in Jackson to derive an explicit   formula for the electric field due to a dipole. Show    numerically, that for large distances (take a couple    of random observation points) the potential function    derived analytically, gives nearly the same values as  those calculated from the exact expression.";