F:=y^3-3*x^2*y^2+3*x^4*y-x^6-x^7*y+x^9+0.000001;fsolve(discrim(F,y),x,complex);fsolve(subs(x=0,F),y,complex);epx:=0.5:epy:=0.8:assume(t,real):assume(u,real):assume(v,real):cube:=subs(x=t+epx*I,y=u+v*I,F):Fr:=Re(expand(cube)):Fi:=Im(expand(cube)):Iv:=subs(x=-epx+I*epx,F):Lo:=fsolve(Iv,y,complex);Now we need to select those with modulus less then epy.Iv1u:=Re(Lo[1]):Iv1v:=Im(Lo[1]): Iv2u:=Re(Lo[2]):Iv2v:=Im(Lo[2]): Iv3u:=Re(Lo[3]):Iv3v:=Im(Lo[3]):Tr:=diff(Fr,t)+diff(Fr,u)*uprime+diff(Fr,v)*vprime: Ti:=diff(Fi,t)+diff(Fi,u)*uprime+diff(Fi,v)*vprime:ode:=subs(u=U(t),v=V(t),uprime=diff(U(t),t),vprime=diff(V(t),t),Tr), subs(u=U(t),v=V(t),uprime=diff(U(t),t),vprime=diff(V(t),t),Ti):fncs:={U(t),V(t)}:cur1:= dsolve({ode,U(-epx)=Iv1u,V(-epx)=Iv1v},fncs,type=numeric,method=classical): cur2:= dsolve({ode,U(-epx)=Iv2u,V(-epx)=Iv2v},fncs,type=numeric,method=classical): cur3:= dsolve({ode,U(-epx)=Iv3u,V(-epx)=Iv3v},fncs,type=numeric,method=classical):with(plots):frm:=pointplot3d({[epx,epy,epy],[epx,epy,-epy],[epx,-epy,epy],[epx,-epy,-epy], [-epx,epy,epy],[-epx,epy,-epy],[-epx,-epy,epy],[-epx,-epy,-epy]},color=black):pl1:=odeplot(cur1, [t,U(t),V(t)], -epx .. epx, numpoints=100,thickness=3,color=red): pl2:=odeplot(cur2, [t,U(t),V(t)], -epx .. epx, numpoints=100,thickness=3,color=blue): pl3:=odeplot(cur3, [t,U(t),V(t)], -epx .. epx, numpoints=100,thickness=3,color=green):display(frm,pl1,pl2,pl3,axes=boxed,scaling=constrained,tickmarks=[0,0,0]);