Consider doing what I did when you first posted this.  That is, wrap your code into a Maple procedure, then use the Maple debugger to step through the code.  That is, do

myproc := proc()
  (* all your code *)
end proc:

stopat(myproc):
myproc();

That will bring up the Maple debugger.  You can execute each line one at a time, using the next, into, step, etc commands (read the ?debugger help page).

The problem is equivalent to partitioning a set of 2*m objects into m-partitions of 2 objects.  That can be easily done recursively:

part2 := proc(L :: list)
local n,k,p;
    n := nops(L);
    if n = 2 then
        return [[L]];
    end if;
    return [seq(seq([ [L[1], L[k]], p[] ]
                     , p in procname(subsop(1=NULL,k=NULL,L)))
                , k=2..n )]
end proc:

 part2([seq(1..4)]);
            [[[1, 2], [3, 4]], [[1, 3], [2, 4]], [[1, 4], [2, 3]]]

The conversion of the partitioned integers to the desired form is straightforward. 

Another possibility, which allows only polynomials, and not, say, a+sin(b), is

select(type, L, polynom(constant, {a,b,c}));

That will allow a + b as well as Pi, which may or may not be desired.

See ?type,polynom.

Hmm.  Where does t enter the expressions?  Using all numeric values my machine executes that nested loop (P=100, Q=20, R=20) in about 0.75 seconds.  Using evalhf (and declaring the arrays appropriately) it does it in 0.2 seconds.  Leaving all non-arrays as symbolic it took 2.2 seconds.  This is using an i5 processor.

Note that you can use ?seq with the optional step argument to specify the rows:

 A := Matrix(10,2, (i,j)-> cat(a,i,j));          
                                   [a11     a12 ]
                                   [a21     a22 ]
                                   [a31     a32 ]
                                   [a41     a42 ]
                             A :=  [a51     a52 ]
                                   [a61     a62 ]
                                   [a71     a72 ]
                                   [a81     a82 ]
                                   [a91     a92 ]
                                   [a101    a102]

A[[seq(1..10,2)], ..];                
                                 [a11    a12]
                                 [a31    a32]
                                 [a51    a52]
                                 [a71    a72]
                                 [a91    a92]

No, there isn't a good way to do so.  This sounds like a bug.  You might be able to work around it by saving the model, exiting, then restarting MapleSim.  Let me know if that works.

Here's one approach, using ?selectremove to partition the solutions into two lists.  Note that I used ?fsolve rather than solve; that seems appropriate considering the equation contains floating-points (i.e. is not exact).

det := eval(Determinant(eigen), lambda=1);
fsol := [fsolve(det, m, 'complex')];
(re,cx) := selectremove(type, fsol, realcons);

The immediate problem occurs during the first loop:

x4a, x4b := fsolve((x-h)^2+(eq4-k)^2-r^2);

The argument to fsolve evaluates to

(x-2)^2+87.04

which has no real solutions, hence the assignment error.  You could specify the complex option to ?fsolve, but do you want complex solutions?

This is discussed in the ?algsubs help page; see the last paragraph of the second bullet.  You can do

 foldr(algsubs, q_x, rho=rho_an, 1/rho=1/rho_an);
                         / /d   \           /d       \\
                       k |-|-- p| alpha + p |-- alpha||
                         \ \dx  /           \dx      //
                       --------------------------------
                                          2
                                   R alpha

Try ?PolynomialTools[CoefficientVector].

Here's a module that you can use to analyze the Parrondo paradox.

Parrondo := module()
option package;
export ExpectedValue
    ,  Pstationary
    ,  ProbOfEvents
    ,  A, B;
    ;
global e;

    # Define the state-transition matrices.
    # Pij is the probability of going from
    # state i to state j.  
    
    A := Matrix(3, [[   0   , 1/2-e , 1/2+e ],
                    [ 1/2+e ,   0   , 1/2-e ],
                    [ 1/2-e , 1/2+e ,   0   ]]
               );
    
    B := Matrix(3, [[   0   , 1/10-e , 9/10+e ],
                    [ 1/4+e ,   0    , 3/4-e  ],
                    [ 3/4-e , 1/4+e  ,   0    ]]
               );
    
    Pstationary := proc( P :: list(Matrix) )
    local lambda,M,V,y,Y,pos;
    description "Compute the stationary probability vector";
        # M is the transition matrix of the system.
        M := foldl(`.`, P[]);
        # The stationary probability vector is a row vector 
        # that satisfies Ps.M = Ps.  Use Eigenvectors
        # with the transpose of M to compute the tranpose of Ps.
        (lambda,V) := LinearAlgebra:-Eigenvectors(M^%T);
        if not member(1, convert(lambda,list), 'pos') then
            error "cannot find unity eigenvalue";
        end if;
        Y := V[..,pos];
        # Normalize the vector.
        return normal~(Y/add(y, y in Y));
    end proc;
    
    ExpectedValue := proc( P :: list(Matrix) )
    local i,k,bits,evts,val,E,n,Ps;
    description "Compute the expected value";
        Ps := Pstationary(P);
        n := nops(P);
        E := 0;
        for k from 0 to 2^n-1 do
            # Sum value of each possible sequence of evts.
            bits := Bits:-Split(k,':-bits'=n);
            evts := subs(0=-1, bits);
            val := add(e, e in evts);
            if val <> 0 then
                E := E + val*add(Ps[i]*ProbOfEvents(i, evts, P), i = 1..3)
            end if;
        end do;
        return normal(E);
    end proc;
    
    ProbOfEvents := proc(i0::{1,2,3}, evts::list({-1,+1}), P :: list(Matrix))
    description "Compute the probability of a given sequence of events, starting at state i0";
    local i,j,k,p;
        i := i0;
        p := 1;
        for k to nops(evts) do
            j := i + evts[k];
            j := 1 + modp(j-1, 3);
            p := p*P[k][i,j];
            i := j;
        end do;
        return p;
    end proc;
    
end module:

For example, to compute the expected (stationary) value of the game A, A, B, B, do

with(Parrondo):
E := ExpectedValue([A,A,B,B]);
eval(E, e=0);
                                           16/163
plot(E, e=0..0.02);

Be careful with has(M,1).  It is not documented to work only on the entries of rtables, though it appears to.  That is, from the help page one might conclude that

has(Matrix(1,1,[[3]]) 

should return true because one of its subexpressions from ?op does:

op(Matrix(1,[3]));        
                1, 1, {(1, 1) = 3}, datatype = anything, storage = rectangular, order = Fortran_order, shape = []

That is not the case, but that doesn't necessarily mean it does what you want.  As an example (possibly not relevant):

has(1/x, -1);
                                   true

Earlier this week I wrote sorting with a permutation list, which shows how to do just that. That is, it generates a permutation that can then be applied to x to get y. Applying it here

x := [1,5,3,7]:
P := map(attributes, sort([seq(setattribute(evalf(x[i]),i), i = 1..nops(x))]));
                         P := [1, 3, 2, 4]

[seq(x[P[i]], i = 1..4)];
                         [1, 3, 5, 7]

You could avoid the issue by using rationals, say 1/10 or 1/20.

If Maple would display an overloaded operator in infix rather than prefix form, then it would be easy:

delayadd := module()
option package;
export `+`;
end module:

with(delayadd):
1+2+3+4;
                           +(+(+(1, 2), 3), 4)

You could do

delayadd := module()
option package;
export `+`;
    `+` := proc() :-`+`(map(``, [_passed])[]) end proc;
end module:

with(delayadd):
1+2+3+4;
                         ( ( (1) +  (2)) +  (3)) +  (4)

but the resulting nesting isn't very nice