coupled second order differential equation

edit

If you want a numerical solution, you can try the following.

a=0.5
var('x,y,u,v,t')
tspace=srange(0,10,.05)
dx1 = u
du = -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*v
dy1 = v
dv = -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*u
result=desolve_odeint([dx1,du,dy1,dv],[1,0,1,0],tspace,[x,u,y,v],ivar=t)
xy=zip(result[:,0],result[:,2])
tx=zip(tspace,result[:,0])
ty=zip(tspace,result[:,2])
tu=zip(tspace,result[:,1])
tv=zip(tspace,result[:,3])

html('t vs. x,y')
show(line(tx)+line(ty,color='green'),axes_labels=['$t$','$x,y$'])
html('phase plane (x vs. y)')
show(line(xy),axes_labels=['$x$','$y$'],aspect_ratio=1)
html('t vs. u,v')
show(line(tu)+line(tv,color='green'),axes_labels=['$t$','$u,v$'])


fu(x,u,y,v) = -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*v
fv(x,u,y,v) = -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*u


tudot=[(tspace[i],fu(result[i][0],result[i][1],result[i][2],result[i][3])) for i in range(0,len(tspace))]

tvdot=[(tspace[i],fv(result[i][0],result[i][1],result[i][2],result[i][3])) for i in range(0,len(tspace))]

show(line(tudot,color='red')+line(tvdot,color='orange'),axes_labels=['$t$','$\dot{u},\dot{v}$'])

You can run the code at the single cell server by clicking here.

You'll need to let Sage know what the independent variable is by including ivar=t in your command. Further, since you are solving a system of equations, you need to use desolve_system.

a,t = var('a,t')
x = function('x', t)
y = function('y', t)
dx2 = 0 == -diff(x,t,2)-(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*diff(y,t,1)
dy2 = 0 == -diff(y,t,2)-(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*diff(x,t,1)
desolve_system([dx2,dy2],[x,y],ivar=t)

This code gets farther but gives an error as well. This might be related to limitations in Maxima's inverse Laplace transform abilities.