Revision history [back]

I don't think there's a way to automate this totally. 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)

I don't think there's a way to automate this totally. 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)