1 | initial version |

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))
```

2 | changes to reflect edit to original question |

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))

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.