with(Statistics):
with(plots):
#Same incoming position but energy first negative then positive
G:=1: x[0]:=-1: y[0]:=1: dx[0]:=.8: dy[0]:=0:  h:=.1: m:=115: 
#G:=1: x[0]:=-1: y[0]:=1: dx[0]:=1: dy[0]:=0:  h:=.1: m:=300: 
#G:=1: x[0]:=-1: y[0]:=1: dx[0]:=1.1: dy[0]:=0:  h:=.1: m:=700: 
#G:=1: x[0]:=-1: y[0]:=1: dx[0]:=1.19: dy[0]:=0:  h:=.1: m:=100: 
#G:=1: x[0]:=-1: y[0]:=1: dx[0]:=1.25: dy[0]:=0:  h:=.1: m:=190: 
#G:=1: x[0]:=-1: y[0]:=1: dx[0]:=1.3: dy[0]:=0:  h:=.1: m:=100: 
E[0]:=evalf((1/2)*(dx[0]^2+dy[0]^2)-G/sqrt(x[0]^2+y[0]^2)):
print("Energy",E[0]):
#print("E",E[0],"gM",G,"position",[x[0],y[0]],"velocity",[dx[0],dy[0]]):
#x[0]:=0: y[0]:=1: dx[0]:=0: dy[0]:=1:  #Rocket
g[1]:=dx:                  #dx/dt=dx
g[2]:=-G*x/((x^2+y^2)^(3/2)):    #d^2x/dt^2=d(dx)/dt
g[3]:=dy:                  #dy/dt=dy
g[4]:=-G*y/((x^2+y^2)^(3/2)):    #d^2y/dt^2=d(dy)/dt
interface(quiet=true):
stoplim:=m:
for n from 0 to m-1 do: 
  for i from 1 to 4 do
    k[n][1][i]:=evalf(h*subs(x=x[n],dx=dx[n],y=y[n],dy=dy[n],g[i])):
  od:
  for i from 1 to 4 do
    k[n][2][i]:=evalf(h*subs(x=x[n]+k[n][1][1]/2,dx=dx[n]+k[n][1][2]/2,
                       y=y[n]+k[n][1][3]/2,dy=dy[n]+k[n][1][4]/2,g[i])):
  od:
  for i from 1 to 4 do
    k[n][3][i]:=evalf(h*subs(x=x[n]+k[n][2][1]/2,dx=dx[n]+k[n][2][2]/2,
                        y=y[n]+k[n][2][3]/2,dy=dy[n]+k[n][2][4]/2,g[i])):
  od:
  for i from 1 to 4 do
    k[n][4][i]:=evalf(h*subs(x=x[n]+k[n][3][1],dx=dx[n]+k[n][3][2],
                        y=y[n]+k[n][3][3],dy=dy[n]+k[n][3][4],g[i])):
  od:

  x[n+1]:=evalf(x[n]+(1/6)*(k[n][1][1]+2*k[n][2][1]+2*k[n][3][1]+k[n][4][1])):
  dx[n+1]:=evalf(dx[n]+(1/6)*(k[n][1][2]+2*k[n][2][2]+2*k[n][3][2]+k[n][4][2])):
  y[n+1]:=evalf(y[n]+(1/6)*(k[n][1][3]+2*k[n][2][3]+2*k[n][3][3]+k[n][4][3])):
  dy[n+1]:=evalf(dy[n]+(1/6)*(k[n][1][4]+2*k[n][2][4]+2*k[n][3][4]+k[n][4][4])):
  E[n+1]:=evalf((1/2)*(dx[n+1]^2+dy[n+1]^2)-G/sqrt(x[n+1]^2+y[n+1]^2)):
od:

X:=[seq([x[n],y[n]],n=0..stoplim)]:
clis:=[]:
for i from 1 to stoplim do clis:=[op(clis),blue] od:
plot(X,color=red):
pointplot(X,color=[red,op(clis)]);