How to plot field lines in sage

edit

I don't think there's a way to automate this totally and avoid putting in all of the initial conditions. Here is another version.

var('x y')
E1=vector([x,3*y])/(sqrt(x^2+y^2))^(3/2)
E2=E1.subs(x=x-1)
F=(E1-E2)/norm(E1-E2)
int_conds=[[0.05,0.001],[0.05,0.005],[0.05,-0.001],[0.05,-0.005]]
solns=[desolve_odeint(F,init,srange(0,1.5,0.05),[x,y]) for init in int_conds]
nhalf=len(solns[0])/3
a=[arrow(solns[i][nhalf],solns[i][nhalf+1]) for i in range(0,len(solns))]
b=map(lambda x: line(x,thickness=2),solns)
show(sum(a)+sum(b)+point((0,0),size=50)+point((1,0),size=50)+plot_vector_field(F,(x,-0.1,1.1),(y,-.5,.5)),axes=false)

Thanks for the edit. See if something like this does what you want.

var('x y')
E1=vector([x,y])/(sqrt(x^2+y^2))^3
E2=E1.subs(x=x-1)
F=(E1-E2)

tstop=0.4
step=0.005
int_conds=[[0.5,0.1*k] for k in [1..3,-3..-1]]
solns=sum([line(desolve_odeint(F,init,srange(0,tstop,step),[x,y]),thickness=2) for init in int_conds])
solns2=sum([line(desolve_odeint(F,init,srange(0,-1*tstop,-1*step),[x,y]),thickness=2) for init in int_conds])

F=(E1-E2)/norm(E1-E2)
solns+solns2+plot_vector_field(F,(x,-0.2,1.2),(y,-1,1))

In the attached file one can see the maxima version (code in my previous answer)

sdfplot.pdf

I'm not sure whether you already have equations to use or not. How you plot these all depends, of course, on how the equations are given. But, here is something. Using information in section D.73.3 of this webpage on quantum mechanics, we need to plot ${(x^2+z^2)}^{(3/2)}=c x^2$ for various values of $c$. The following code produces a plot.

var('x z')
sum([implicit_plot((x^2+z^2)^(3/2)==c*x^2,(x,-5,5),(z,-5,5)) for c in [1,2,3,4]])

People interested in this question will want to perhaps help out in making our ticket for streamline plots a reality.