@tomleslie 

The problem with all three of your methods is that they cause the limit to be evaluated multiple times. You can give the command debug(limit) before executing G to see how many times limit is called. In addition to just being generally wasteful, it is clear in the Question (and even in its title) that the OP wants the limit evaluated before the call to G (indeed, even before the definition of G).

You may be especially surprised that your third method also suffers from multiple limit evaluations; however, it does because of the special evaluation rules of add (which also apply to seq and mul), which re-evaluate the first argument for every value of the index. So, you may think that you can avoid this by calling limit before the add but still inside G. Something like this:

G:= proc(ff::procedure, x, y, n)
local ll, L:= limit(ff(x,y), y= ll);
     add(L, ll= 1..n)
end proc:

This does call limit only once, but the answer returned is wrong, indeed, ridiculous, as it still contains ll. I urge you to try it yourself and see.

The way around this is to use subs to modify a "shell" procedure into a working procedure, like this:

f:= (x,y)-> sin((x+y)/y):
F:= (f,x,y,L)-> limit(f(x,y), y= L):
G:= proc(x, y, n)
local ll;
     add(_FF, ll= 1..n)
end proc:

G1:= subs(_FF= F(f,x,y,ll), eval(G)):

G1(a,b,10);

Note that this allows reference to the locals and parameters of G outside of G itself. 

There is another "clever" way to achieve this with unapply and no subs:

G:= unapply('`+`'(F(f,x,y,ll) $ ll= 1..n), [x,y,n]);

Note that the + is enclosed in both name quotes and unevaluation quotes.

@w-Teilchen Use solve---it's even more direct than roots, as in -solve(numer(q)), etc.

Forget about patmatch. It's definitely not the way to solve your problem.

If both solve and roots don't work for you, please post an example problem for which that is true.

@wolverine Okay, you've substantially cleaned up the code by using my redefinition of ode. But you haven't answered my question about what you mean by D[1] and D[2]. As it stands, they are being treated as the derivatives with respect to the first and second independent variables, respectively. Is that what you intend?

@MDD Please post this as a new Question (a new thread). This is far, far removed from the original Wronskian-related Question of this thread. It has very little to do with PolyLinearCombo. It's substantially more complicated than your previous questions. And I think that more people will seriously think about it if it's a new Question.

In the new Question, please mention

1. whether f and the members of F are always homogenous in the variables. (I don't immediately see that that makes a difference, but any restriction that you can make on the complexity of the problem can't hurt.)
2. whether the coefficients of f, the coefficients of the members of F, and the members of N are all always polynomials in the parameters (or contant). 
   

@Dmitry Lyakhov Mac Dude's nested loop code should work. If it doesn't, perhaps you have an installation bug. Please post an executed worksheet showing it not working.

@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;