@acer Many thanks.  Would it be possible for you to provide a small example of the implementation you have in mind, particularly the approach based on tables (or another more suitable data structure)? I would greatly appreciate seeing a concrete implementation, as this would help me better understand how to adapt the procedure for large-scale computations. I do apologize for asking this again and for not having learned from the countless times you and Carl Love have helped me in the past.

@janhardo Thank you.   

@janhardo Pretty impressive. Thank you very much. Unfortunately, I tend to overcomplicate things when writing a procedure for the first time. Also, Vars is no longer needed.

@C_R There is no problem with the uploaded example itself. However, when I run the procedure over a large pool of models (for example, systems with 6 monomials instead of 5, resulting in about 203685 systems in total), Maple eventually stops progressing at some point. At that stage, CPU usage drops from nearly 100% across all cores to almost zero, and Maple appears to stall indefinitely without producing an error message.

I came up with the following simplified procedure that passes the ThreadSafetyCheck,

generateNonlinearModelsPlus_fast := proc(model::list, fullmodel::list, vars::list := [x,y,z])
description "Generate all models with one additional monomial from fullmodel":
local n, i, j, diff, coeff, mon, out, k;

    n := nops(model);
    out := table();
    k := 0;

    for i to n do
        # One expansion is enough; simplify is usually much more expensive here.
        diff := expand(fullmodel[i] - model[i]);

        if diff <> 0 then
            coeff := [coeffs(diff, vars, 'mon')];
            mon := [mon];

            for j to nops(coeff) do
                k := k + 1;
                out[k] := subsop(i = model[i] + coeff[j]*mon[j], model);
            end do;
        end if;
    end do;

    return [seq(out[i], i = 1 .. k)];
end proc:

@sand15 

The idea is to generate all systems containing a total of (n+1) monomials from a pool of systems containing (n) monomials that satisfy a given condition. Each new system is obtained by adding a single monomial selected from the set difference between the full system and the corresponding (n)-monomial system.

As an example, the full system may correspond to a 3D polynomial system containing all possible monomials up to degree 3.

@sand15 Many thanks. Here it is an example of how to run the procedure.

Example_Generate_Models.mw

As for threadsafety, I knew that zip was not threadsafe, thus using Grid.

@acer Thank you for the explanation. fixup saves me a lot of time when transferring my results from Maple to LaTeX.  

@acer Thank you, but I think I may have missed something. fixup returns the same input expression expr, doesn't it? How was it modified?

@acer Sorry to reopen the thread, but do you happen to know how to fix z*x? I mean…

from

[sigma*y, (-sigma - 1)*y + z*x, (rho - 1 - z)*beta - (x + y)*x]

to

[sigma*y, (-sigma - 1)*y + x*z, (rho - 1 - z)*beta - (x + y)*x]

@acer I could not ask for more.  Thank you.  

@acer I could never have come up with such a clean procedure—it's got Carl written all over it. I can't thank him enough for all the help he's given me. And now, you're part of that too.

Example

lorenz := [-sigma*x + sigma*y, rho*x - x*z - y, -beta*z + x*y];
vars := [x,y,z];

map(w->monomialsGB(w,vars),lorenz);

and output is

[[sigma, -sigma], [y, x], [-1, -1, rho], [z*x, y, x], [-beta, 1], [z, y*x]]

The ordering should be [x,y], [x,y,x*z] and [z,x*y].  

Actually I need both versions. I mean Carl's and yours.   

Many thanks.

@acer Thank you! That’s perfect. I came across only a few cases where repetition occurs, and this will save me the trouble of handling them separately.

Admitting my total lack of skill (even though I've been a Maple user for years) and hoping it's not too much trouble — how could the following procedure be modified to use the ordering established above?

monomialsGB := proc (p::algebraic, vars::list)
local
     M,
     v:= indets(p, suffixed({vars[]})),
     C:= coeffs(expand(p), v, 'M'),
     P:= plex(ListTools:-Reverse([v[]])[]),
     S:= sort(`[]`~([C],[M]), (a,b)-> not Groebner:-TestOrder(a[2], b[2], P));
     (map2(op, 1, S), map2(op, 2, S)):
end proc:

Please include details on how you modified the procedure. I really need to learn the basics of Maple once and for all.

@acer Many thanks.  Yes, there is the possibility of repeated terms.

@acer As usual, you went above and beyond—thank you for that! It's a great way for me to learn more about Maple commands. In the cases I typically work with, the entries usually have the same level of nesting.

@acer Thank you for making it clear.  The inner most entries are number (integers, rationals and possbily reals) and/or symbols (variables), or empty. In cases I am interested in they are a solution of an equation.