EOM[x0_,y0_,p0_,q0_]:={x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]y[t]^2-A x[t],q'[t]==-y[t]x[t]^2-A y[t],x[0]==x0,y[0]==y0,p[0]==p0,q[0]==q0}

H[x_,y_,p_,q_]:=p^2+q^2+x^2 y^2+A(x^2+y^2)

qpos[x_,y_,p_]:=Sqrt[energy-p^2-x^2 y^2-A(x^2+y^2)]

flow[{xi_,yi_,pi_,qi_},T_]:=
{x[T],y[T],p[T]}/.NDSolve[EOM[xi,yi,pi,qi],{x,y,p,q},{t,0,T}][[1]]

render[g_]:=Map[Polygon,Join[Flatten[Table[{g[[i,j]],g[[i,j+1]],g[[i+1,j]]},{i,n},{j,n}],1],Flatten[Table[{g[[i+1,j+1]],g[[i,j+1]],g[[i+1,j]]},{i,n},{j,n}],1]]]

time=500;
A=.25;
x0=1;
y0=0;
p0=0;
q0=1;
energy=H[x0,y0,p0,q0]

trajectory=NDSolve[EOM[x0,y0,p0,q0],{x,y,p,q},{t,0,time}];

orbit=ParametricPlot[Evaluate[{x[t],y[t]}/.trajectory],{t,0,time},PlotRange->All,Axes->False];

contour=ContourPlot[H[x,y,0,0]==energy,{x,-2.5,2.5},{y,-2.5,2.5},ContourStyle->Red];

Show[orbit,contour]

w=.1;
n=6;
f1=Table[{x0,y0+w(i/n),p0+w(j/n),qpos[x0,y0+w(i/n),p0+w(j/n)]},{i,0,n},{j,0,n}];
f2=Table[{x0+w,y0+w(i/n),p0+w(j/n),qpos[x0+w,y0+w(i/n),p0+w(j/n)]},{i,0,n},{j,0,n}];
f3=Table[{x0+w(i/n),y0,p0+w(j/n),qpos[x0+w(i/n),y0,p0+w(j/n)]},{i,0,n},{j,0,n}];
f4=Table[{x0+w(i/n),y0+w,p0+w(j/n),qpos[x0+w(i/n),y0+w,p0+w(j/n)]},{i,0,n},{j,0,n}];
f5=Table[{x0+w(i/n),y0+w(j/n),p0,qpos[x0+w(i/n),y0+w(j/n),p0]},{i,0,n},{j,0,n}];
f6=Table[{x0+w(i/n),y0+w(j/n),p0+w,qpos[x0+w(i/n),y0+w(j/n),p0+w]},{i,0,n},{j,0,n}];

T=.1;
g1=Map[flow[#,T]&,f1,{2}];
g2=Map[flow[#,T]&,f2,{2}];
g3=Map[flow[#,T]&,f3,{2}];
g4=Map[flow[#,T]&,f4,{2}];
g5=Map[flow[#,T]&,f5,{2}];
g6=Map[flow[#,T]&,f6,{2}];

plot=Graphics3D[{EdgeForm[Black],Map[render,{g1,g2,g3,g4,g5,g6}]},Boxed->False]

(* T=6 *)

(* T=8 *)

(* T=10 *)

(* T=12 *)

(* T=14 *)