
V=1/Sqrt[(x-1)^2+y^2]+1/Sqrt[(x+1)^2+y^2]
m = 5
            
Plot3D[V,{x,-m,m},{y,-m,m},PlotRange->{{-m,m},{-m,m},{0,10}}]

Ex=-D[V,x]
Ey=-D[V,y]
            
Plot[Ey/.{x->0},{y,-m,m},AxesLabel->{"y","Ey"}]
StreamPlot[{Ex,Ey},{x,-m,m},{y,-m,m},StreamPoints->Fine,AxesLabel->{"x","y"}]

b=-1
x0=0.2
y0=5
m=5
equation={x'[t]==vx[t],vx'[t]==b*Ex/.{x->x[t],y->y[t]},y'[t]==vy[t],vy'[t]==b*Ey/.{x->x[t],y->y[t]},x[0]==x0,vx[0]==0,y[0]==y0,vy[0]==0}
tmax = 50
sol = NDSolve[equation,{x[t],vx[t],y[t],vy[t]},{t,0,tmax}]
            
ParametricPlot[{x[t],y[t]}/.sol,{t,0,tmax},PlotPoints->1000,PlotRange->{{-m,m},{-m,m}},AspectRatio->1]

tmax = 200
sol = NDSolve[equation,{x[t],vx[t],y[t],vy[t]},{t,0,tmax}]
             
ParametricPlot[{x[t],y[t]}/.sol,{t,0,tmax},PlotPoints->1000,PlotRange->{{-m,m},{-m,m}},AspectRatio->1]

tmax=4000
sol = NDSolve[equation,{x[t],vx[t],y[t],vy[t]},{t,0,tmax},
                Method->{"ExplicitRungeKutta","DifferenceOrder"->5},
                AccuracyGoal->6,PrecisionGoal->6]
              
Plot[x[t]/.sol,{t,0,tmax},PlotPoints->1000,AxesLabel->{"t","x"},PlotRange->{{0,tmax},{-5,5}}]

energie=1/2*(vx[t]^2+vy[t]^2)-V/.{x->x[t],y->y[t]}
              
Plot[energie/.sol,{t,0,tmax},PlotPoints->1000,AxesLabel->{"t","x"},PlotRange->{{0,tmax},{-0.5,0}}]
