@vv It's feasible to do the 2x2 case for moduli q up to about 64.

I'd forgotten that Iterator was now available in out-of-the-box Maple. It's great that there's finally a library alternative to combinat:-cartprod.

As far as I know, RootFinding:-Isolate is the only sure way to get (approximations of) all of the solutions for a polynomial system.

@Kitonum You've come up with a brilliant new way to do Cartesian products in Maple! There hasn't been a new good one (that I've seen) in nine years. You should add this procedure as a comment to this 2007 thread: Cartesian product of lists. Maybe Joe Riel will incorporate it into the post. The parameter Expr should be typed uneval instead of algebraic. Then you won't need the unevaluation quotes, and your NestedSeq will be more like regular seq.

However, your code above is about 50% slower than my most recently posted for loop. I think that in general a for loop is faster than a seq when then seq has a significant percentage of NULLs. Still, your procedure will be great for cases where there's not a significant percentage of NULLs.

@vv Yes, using LinearAlgebra:-Determinant is faster, and LinearAlgebra:-Modular:-Determinant is even faster than that (by a little bit). Here's the modified code:

restart:
macro(LA_M= LinearAlgebra:-Modular):
Wordsize:= integer[kernelopts('wordsize')/8]:
CartProdSeq:= proc(L::seq(list), $)
local S,j;
    eval(subs(S= seq, foldl(S, [_i||(1..nargs)], seq(_i||j= L[j], j= 1..nargs))))
end proc:

PRIME:= 2:  n:= 4:

for M in CartProdSeq([CartProdSeq([$0..PRIME-1] $ n)] $ n) do
     if LA_M:-Determinant(PRIME, LA_M:-Mod(PRIME, Matrix(M), Wordsize)) <> 0 then
          GROUP[cat("", op~(M)[])]:= [][]
     end if
end do:

GROUP:= [indices(GROUP, nolist)]:

@acer Yes, that's my current setup almost exactly: Maple 16.02, 64-bit, Windows 7, Standard GUI.

@Stavros Sorry, I needed to use Det instead of Determinant. No packages are needed. Here's the corrected code:

restart:

CartProdSeq:= proc(L::seq(list), $)
local S,j;
    eval(subs(S= seq, foldl(S, [_i||(1..nargs)], seq(_i||j= L[j], j= 1..nargs))))
end proc:

PRIME:= 2:

for M in CartProdSeq([$0..PRIME-1] $ 16) do
     if Det(Matrix((4,4), M)) mod PRIME <> 0 then GROUP[cat("", M[])]:= [][] fi
end do:

GROUP:= [indices(GROUP, nolist)]:
nops(GROUP);

     20160

@acer That's great, especially the resolution. What does the color indicate? Is it the number of attractors per r value?

@Joe Riel I can't reproduce it either anymore. But I reproduced it about a dozen times before posting.

However, there's still the issue of the args being reported as a large integer by trace. That's not a big concern, but I wonder if it's happening for you.

I just finished the code and posted it to the "Bifurcation in Maple" thread with some very interesting plots.

@Joe Riel There are often requests here for some way to circumvent automatic simplification for display purposes. Your 3/3 shows that you have done this. I've also checked that it is possible to make the denominator 0. These forms seem to be stable. I got 3/3 + 3 = 12/3, and 1/0 + 1 = 1/0.

I have often wished that evalindets had an option recurse which would cause it to use eval[recurse] instead of regular eval.

How do you use ifelse? I can find no documentation, nor does eval(ifelse) return anything useful.

@Les When you're entering or editing the Question, there's a pull-down menu that lists Maple, MapleSim, MapleTA, MaplePrimes, etc.---the various Maplesoft products that are discussed on this forum. While the vast majority of the discussion here is about Maple itself, there are occasionally Questions about the others.

@Les The single quotes are used to guard against the possibility that elsewhere in your code you assigned a value to color, linestyle, spacedash, etc. If you've not assigned values to any of those, then the quotes make no difference. If you want to get even more fussy about it, you can make them, for example, ':-color', which would also guard against the possibility that the code is embedded in other code where color is used as a local variable.

To deal with indexed names and functions in the most general way, I think that it's necessary to localize the names (symbols) before they are assigned values. Like this:

RootSymbol:= (e::{symbol, indexed, function})-> `if`(e::symbol, e, thisproc(op(0,e))):
     
Localize:= proc(e)
local S:= RootSymbol~(indets(e, {symbol, indexed, function}));
     subs(S =~ convert~(S, `local`), e)
end proc:

Define your list or vector of column headers/entries:
Exprs:= [a, b(23)[34], c[b[a]^2+1]];

Headers:= Localize(Exprs):
a:= 1:  b:= B:  c:= f:
Matrix([Headers, Exprs]);

Now you can change the values of a, b, c at will, and the Headers will never change.

@rostam And do you have the same anti-aliasing setting on your laptop and PC? Use the menu Tools -> Options -> Display -> Plot anti-aliasing.