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

You can do substantially better, but with a bit of work.  First, let's compute an expression for the number of valid lists.  Given the integers 1..n, we know that the minimum allowable k (number of zeros) is n-1.  Let z be the number of excess zeros, that is, z = k-n+1.  For the integers 1..n-1, each must be followed by a zero.  So think of the integers as the pairs (1,0), (2,0), ... (n-1,0) and the singleton n.  We have z zeros remaining to place among these n objects.  This is equivalent to choosing z objects from n+z objects, or binomial(n+z,z) = binomial(k+1, k-n+1).  So for the n=6, k=8 example, the number of lists is binomial(9,3) = 84.

To make this sort of thing really fast, in Maple, consider compiling the code.  Furthermore, to reduce the memory usage and gain speed, generate and operate on each output list one at at a time, rather than generating all the lists at once.  More accurately, rather than generating a list of lists, replace, in a loop, the contents of a mutable data structure (an rtable).

I have written a module that I hope to release soon (still needs some cleanup) that creates fast (compilable) iterators for common finite structures. I'll demonstrate it here. The procedure that generates each of the structures you want is

myp1 := proc(n,k)
local z, P,A,Z,W;
    z := k-n+1; # excess zeros
    (P,A,Z) := Iterator:-Combination(n+z,z,'compile');
    W := Array(1..n+k, 'datatype' = integer[4]);  # this Array stores the desired result
    while P(A) do  
        LocToComb1(n,z,Z,W);
         printf("%d\n", W);  # normally one would do something useful with the generated Array at this point
    end do;
end proc:

Procedure to convert the Array Z, which contains the logical positions of the excess zeros, to the final form, which is stored in W.

LocToComb1 := proc(n,z, Z :: Array(datatype=integer[4]), W :: Array(datatype=integer[4]))
option autocompile;
local i,j,k,p,posz;
    p := 1;  # index into W
    j := 1;  # index into Z
    k := 1;  # integer index
    posz := Z[j]+1;
    for i to n+z do
        if i = posz then
            # insert excess zero
            W[p] := 0;
            p := p+1;
            if j < z then
                j := j+1;
                posz := Z[j]+1;
            else
                # no more excess zeros; finish up integers
                for k from k to n-1 do
                    W[p] := k;
                    W[p+1] := 0;
                    p := p+2;
                end do;
                W[p] := n;
                break;
            end if;
        else
            # insert integer
            W[p] := k;
            if k < n then
                # inset paired zero
                W[p+1] := 0;
                p := p+2;
                k := k+1;
            else
                # no more integers; finish up
                for p from p+1 to 2*n-1+z do
                    W[p] := 0
                end do;
                break;
            end if;
        end if;
    end do;
    return 0; # not used other than to define return type for compiler
end proc:

Here is an example,
 

myp1(3,4);
0 0 1 0 2 0 3
0 1 0 0 2 0 3
1 0 0 0 2 0 3
0 1 0 2 0 0 3
1 0 0 2 0 0 3
1 0 2 0 0 0 3
0 1 0 2 0 3 0
1 0 0 2 0 3 0
1 0 2 0 0 3 0
1 0 2 0 3 0 0

How fast is it?

time(myp1(6,18)); # with printing disabled
                                       0.140

Contrast this with your procedure

time(myp2(6,18));
                                      5.256

Followup

A simpler way to do this is to generate all permutations of n ones and z zeros, then substitute the k'th 1 with the integer pair (k,0), omitting the zero if at the last position.  This can be done with

PermToPattern := proc(n :: integer, z :: integer
                      , Z :: Array(datatype=integer[4])
                      , W :: Array(datatype=integer[4])
                     )
local i,j,p;
option autocompile;
    j := 1;
    p := 1;
    for i to n+z do
        if Z[i] = 0 then
            W[p] := 0;
            p := p+1
        else
            W[p] := j;
            j := j+1;
            if j > n then
                # done with integers; fill remaining 0's
                for p from p+1 to 2*n+z-1 do
                    W[p] := 0;
                end do;
                return 0;
            end if;
            W[p+1] := 0;
            p := p+2;
        end if;
    end do;
    return 0;
end proc;


myp4 := proc(n,k)
local z, P,A,Z,W;
    z := k-n+1;
    (P,A,Z) := Iterator:-Permute([1$n,0$z],'compile');
    W := Array(1..n+k, 'datatype' = integer[4]);
    while P(A) do
        PermToPattern(n,z,Z,W);
        printf("%d\n", W);
    end do;
end proc:

This is not any faster to run, but it was easier to write the procedure that transforms each permutation into the desired pattern.

time(myp4(6,18));
                                      0.168