Revision history [back]

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/2)
E2=E1.subs(x=x-1)
F=(E1+E2)/norm(E1+E2)
int_conds=[[0.05,0.01],[0.4,0.01],[0.99,0.01]]
solns=sum([line(desolve_odeint(F,init,srange(0,1,0.05),[x,y])) for init in int_conds])
solns+plot_vector_field(F,(x,-2,2),(y,-1,1))

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/2)
E1=vector([x,y])/(sqrt(x^2+y^2))^3
E2=E1.subs(x=x-1)
F=(E1+E2)/norm(E1+E2)
int_conds=[[0.05,0.01],[0.4,0.01],[0.99,0.01]]
solns=sum([line(desolve_odeint(F,init,srange(0,1,0.05),[x,y])) 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])
solns+plot_vector_field(F,(x,-2,2),(y,-1,1))
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))