@patient Inside the loop, include a line such as 

P[k]:= plot(< X | A >):

When the loop is done, to see an animation do

plots:-display(convert(P,list), insequence);

@tomleslie You may have wondered why I didn't mention combinat:-cartprod, especially since I had already used the combinat package in my Answer. The reason is that it's extremely inefficient. See this MaplePrimes post, which compares in great detail several different methods of making a Cartesian product.

@mehran rajabi Y is not a function of X over the full range 0..2*Pi. You need to consider the ranges 0..Pi and Pi..2*Pi separately.

If your goal is simply to plot the parametric curve, then you should simply use a parametric plotting command. There's no need to express Y as a function of X.

plot([X(theta), Y(theta), theta= 0..2*Pi]);

If you want to do it with Y as a function of X, using my procedure YY, then you can merge plots over the two ranges. This can be done within a single plot command:

plot(['YY(XX, 0..Pi)', 'YY(XX, Pi..2*Pi)'], XX= X(0)..X(Pi), color= red);

@MDD You have a good point, and what you suggest was also my first idea. However, not every short monomial order is expressed as a function of its variables, even though all the commonly used ones are. So, as a compromise, I made V the last, optional argument and made its default value be op(T). Here's the updated code:

Homogenize.mw

I guess that the procedure has held up under testing so far?

@Markiyan Hirnyk Very nice. Vote up. I also see that this result doesn't rely on that good luck to which I referred earlier.

@Markiyan Hirnyk 

How does one in practice use these results from eliminate to express Y as a real-valued function of X? (I am not expressing any doubt that it can be done---I have done similar things myself. I just think that it would complete the Answer to see it.)

I believe that it is only good luck (and perhaps lexicographic order of X and Y) that some of these results express theta as a function of X alone.

@Kitonum 

How can it be considered a "solution" when Y is expressed as a function of both X and theta?

@roman_pearce 

I don't see the ambiguity, but perhaps I'm missing something. Certainly your example causes no trouble for my procedure. I ran my procedure on your example in my Answer below.

I don't understand what distinguishes f from the members of F. Isn't f being treated identically to the members of F? So why not just put f in F? If you plan to treat f later as distinguished, just put f in a distinguished position in F, such as first or last, and extract it later.

Your usage of D is highly suspicious. Are you aware that D has a predefined meaning as Maple's differential operator? that D[i] means the derivative with respect to the ith variable of a multivariate function (thus implying that there are two (as yet unnamed) independent variables in your differential equation, making it a PDE, not an ODE)? that you have a variable named Dy where you may have meant D(y), which would make y the dependent variable of the differential equation? that the superscript (2) means not squaring but a double application of the operator? that the y in the first term of your ode1 is being treated as a coefficient of the differential operator rather than its argument?

You should upload your complete code as a downloadable worksheet. However, I can guess what most of the variables are/should be defined as. Let me know if I have this right.

ode:= (D@@2)(y) - epsilon*(1-y^2)*D(y) + y = 0;
eord:= 1;

Good point: Come to think of it, in 16 years of Maple usage, I've never used programmatically the output of the with command. I don't see any use for it other than possibly as a reference for the spelling of the command names, which could just as well be found on the package's help page. Anyway, as suggested by taro, you can use a colon (:) to suppress the displaying of the return value of any command.

@Mac Dude 

In Standard, you can easily make an array plot without using plots:-display. Like this:

P1:= plot(x, x= -2..2):  P2:= plot(x^2, x= -2..2):
_PLOTARRAY(Matrix(2,1, [P1,P2]));

(This only works in the Standard interface. In Classic, plots:-display actually does a significant amount of recomputation to make an array plot. I don't know how this will work with plots automatically sent to files with plotsetup.)

I hope that you can make good use of that.

@Mac Dude 

I've been working on it, but it will be a big job to adapt this to the moderm plot structures and to the Standard GUI. The examples which produce an empty plot in the Standard GUI will produce a correct plot in the Classic GUI. Some examples with fancy tickmarks will give errors in either GUI.

As far as I can see, your expression Branch1Cons doesn't contain any arctan functions, or, indeed, any inverse trigonometric functions at all. So what's the point of all those arctans in the substitution list?

The error I was tracking down was ultimately caused by this plotting difference between Standard and Classic. Execute the following line in both:

plot(x, x= -10..10, view= setattribute(-10., "a banana") .. 10.);

In Classic, it's just the normal expected plot. In Standard, it's an empty pair of axes. This was very hard to track down because of the weird way that negative floats with attributes print (on don't print) in Standard. Try this in both Standard and Classic:

x:= setattribute(-23., "a banana");

This works as far as the kernel is concerned, but it prints nothing in Standard. Even lprint(x) shows nothing.

I decided that this wasn't worth a separate Question, because I doubt that there's any workaround. But if someone finds one, I'd be grateful. This situation breaks the logarithmic plotting package that I have posted (14 years ago) on the Maple Applications Center. That package relies on creating plots with attributed floats in the VIEW.