@sursumCorda Using overload with the seq parameter modifier is a brilliant idea, and I think a quite original idea also. I've expanded on it, in an Answer that I just posted, to get a Maple implementation of Mathematica's Cases (at least to the extent that the OP used in this Question).

@Ronan Disjunctively means that has(e, [x,y]) is equivalent to has(e, x) or has(e, y). If instead it were has(e, x) and has(e, y), then that would be conjunctively.

@charlie_fcl This is the reason that I prefer ifelse (or, equivilently, `if`) over piecewise in this situation: piecewise is a highly symbolic condition-based control structure whose boolean conditions do not need to be evaluated immediately. If the boolean condition can be evaluated immediately, ifelse is much simpler and more efficient.

@Andiguys Here is a completely different way to do the 3D plot. This way does not require the double loop (or any loop):

Cv:= 'Cv': tau0:= 'tau0':
MaxArgs:= op@unapply(
   [TRC(tau2, tau1), `union`(C1,C2,C3)], [tau2, tau1, Cv, tau0]
):
Max:= proc(Cv, tau0)
local tau2, tau1;
    try
        Optimization:-Maximize(
            MaxArgs(tau2, tau1, Cv, tau0), tau2= 0..1, tau1= 0..1,
            assume= nonnegative
        )[1]
    catch: undefined
    end try
end proc
:
plot3d(Max, 50..70, 0.3..0.7, labels= ['Cv', 'tau0', TRC], grid= [5,9]);

The 5 and 9 in the grid option were chosen to correspond to the number and spacing of the values of Cv and tau0 that you used in the double loop. You can increase them if you want.

@Andiguys No, the plot command should not be inside the loop.

In the original double loop, you had s[Cv, tau0]:= Maximize(...). You changed that to s:= Maximize(...). If you change it back to s[Cv, tau0], then you will get the plot.

@charlie_fcl Like this:

(e,f):= ifelse(c<d, [a, 2*b], [3*a, b])[]

You could replace ifelse with piecewise (leaving everything else exactly as is) and get the same result, but IMO ifelse is better.

@charlie_fcl In Maple it can be done as follows, so it's likely to also work in Maple Flow:

(a,b):= (min,max)(c,d)

@Carl Love Okay, I couldn't resist the temptation to make the applyrule work, and I found something that works, although I can't explain why you were getting an infinite loop before.

applyrule(
    x::(anything^And(fraction, 2 &under denom))*y::(Not(1) &under denom)
    = y %* x,
    eq
);

@Carl Love I just added a small change to the code above. I replaced denom(R)=1 with S=1 or denom(R)=1. In other words, if there are no coefficients with denominators OR there are no square roots, then leave the product unchanged.

@SARBAST There is no attached file. Please try uploading it again.

@SARBAST I assume that you want to vary the initial values (x(0), y(0))  over [-6,6] x [-3,6]. We still need initial values z(0) and w(0) in order to do any computation.

@C_R The Wikipedia article "Attractor" contains a definition of basin of attraction. It's the mathematical generalization of the geological / geographical / hydrological concept of drainage basin or watershed.

@C_R I haven't run this yet, but in order to get such a slow and smooth animation to post, the OP very likely included a huge number of extra frames (10 times, or more) in the animation. These extra frames are not needed when viewing in a worksheet because you can slow down the animation with the toolbar controls. Not computing the extra frames should speed up the computation significantly.

@dharr As far as I can tell, the worksheet poly.mw attached to your Answer is the same as when it was originally attached. If this is not true, please include a comment in the worksheet to indicate the new part.

@Andiguys This is to diagnose the problem, not to fix it: Please insert this command after the double loop:

lprint(indices(s, pairs, indexorder));

Then execute the entire worksheet. Then post the entire executed worksheet.