I'll see what I can do.  No promises.  If you do use that code, I recommend changing the Pstationary command to use the NullSpace command, as Robert Israel did in his post. That is, use

    Pstationary := proc( P :: list(Matrix) )
    local M,p,Ps;
    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 NullSpace to
        # solve (M^T-1).Ps^T = 0.
        Ps := LinearAlgebra:-NullSpace(M^%T-1);
        if nops(Ps) <> 1 then
            error "could not find vector";
        end if;
        Ps := Ps[1];
        # Normalize the vector.
        return normal~(Ps/add(p, p in Ps));
    end proc;

Also, consider declaring A and B as globals rather than exports so they can be externally reassigned.  That permits looking for symbolic solutions. For example,

A := Matrix([[  0  ,  a1  , 1-a1 ],
             [1-a2 ,  0   ,  a2  ],
             [ a3  , 1-a3 ,  0   ]
            ]):

B := Matrix([[ 0   ,   b1 , 1-b1 ],
             [1-b2 ,   0  , b2   ],
             [b3   , 1-b3 ,  0   ]
            ]):

Ea := ExpectedValue([A]):
Eb := ExpectedValue([B]):
Eab := ExpectedValue([A,B]):
simplify(Eab, {numer(Ea)=0, numer(Eb)=0});
                                                      0

which indicates that if A and B are fair, than A.B is also fair.  That isn't the case for A.A.B.
 

I'll see what I can do.  No promises.  If you do use that code, I recommend changing the Pstationary command to use the NullSpace command, as Robert Israel did in his post. That is, use

    Pstationary := proc( P :: list(Matrix) )
    local M,p,Ps;
    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 NullSpace to
        # solve (M^T-1).Ps^T = 0.
        Ps := LinearAlgebra:-NullSpace(M^%T-1);
        if nops(Ps) <> 1 then
            error "could not find vector";
        end if;
        Ps := Ps[1];
        # Normalize the vector.
        return normal~(Ps/add(p, p in Ps));
    end proc;

Also, consider declaring A and B as globals rather than exports so they can be externally reassigned.  That permits looking for symbolic solutions. For example,

A := Matrix([[  0  ,  a1  , 1-a1 ],
             [1-a2 ,  0   ,  a2  ],
             [ a3  , 1-a3 ,  0   ]
            ]):

B := Matrix([[ 0   ,   b1 , 1-b1 ],
             [1-b2 ,   0  , b2   ],
             [b3   , 1-b3 ,  0   ]
            ]):

Ea := ExpectedValue([A]):
Eb := ExpectedValue([B]):
Eab := ExpectedValue([A,B]):
simplify(Eab, {numer(Ea)=0, numer(Eb)=0});
                                                      0

which indicates that if A and B are fair, than A.B is also fair.  That isn't the case for A.A.B.
 

Some time ago I recommended the monograph Random Walks and Electric Networks, by Doyle and Snell. It introduces most of these concepts. I don't know whether it is still in print, but see that you can download it

Some time ago I recommended the monograph Random Walks and Electric Networks, by Doyle and Snell. It introduces most of these concepts. I don't know whether it is still in print, but see that you can download it

The list won't contain assignments (a0 := b0); Maple assignments are executed immediately.  Use equations:

eq1 := a0 = b0:
eq2 := a1 = a0+b1:
eq3 := a2 = a1+c0:
eq4 := result = a0+a1+a2:
L := eval([cat(eq, 1..4)]);
                               L := [a0 = b0, a1 = a0 + b1, a2 = a1 + c0, result = a0 + a1 + a2]

 

The list won't contain assignments (a0 := b0); Maple assignments are executed immediately.  Use equations:

eq1 := a0 = b0:
eq2 := a1 = a0+b1:
eq3 := a2 = a1+c0:
eq4 := result = a0+a1+a2:
L := eval([cat(eq, 1..4)]);
                               L := [a0 = b0, a1 = a0 + b1, a2 = a1 + c0, result = a0 + a1 + a2]

 

If you just want to creat the list of equations, without any assignment, then do

[seq(a[i] = b[i]+c[i], i = 1..2)];

If you just want to creat the list of equations, without any assignment, then do

[seq(a[i] = b[i]+c[i], i = 1..2)];

For numeric quantities, Maple compares the numeric values.  For other quantities Maple generally compares addresses. So two lists of numerics are equal only if they are the same list.

For numeric quantities, Maple compares the numeric values.  For other quantities Maple generally compares addresses. So two lists of numerics are equal only if they are the same list.

Interesting kernel bug. The behavior is as if parameter expansion occurred after the evalution.  That is, (in curry)

  subs('_X' = args[2..nargs], () -> p(_X,args))

is being evaluated as

   () -> p(p,args)

then the parameter expansion kicks in and both p's are replaced with f.

While that seems a bug (the expectation is that the parameter expansion is done first), it can be avoided with a modification to curry:

curry := proc(p)
    subs(['_p'=p, '_X'=_rest], () -> _p(_X,args));
end proc:

Interesting kernel bug. The behavior is as if parameter expansion occurred after the evalution.  That is, (in curry)

  subs('_X' = args[2..nargs], () -> p(_X,args))

is being evaluated as

   () -> p(p,args)

then the parameter expansion kicks in and both p's are replaced with f.

While that seems a bug (the expectation is that the parameter expansion is done first), it can be avoided with a modification to curry:

curry := proc(p)
    subs(['_p'=p, '_X'=_rest], () -> _p(_X,args));
end proc:

Admittedly its a toy problem, but for it I'd pass both equations to solve and be done with it:

proc()
local x,y,eq1,eq2;
    eq1 := 5+3*x=0;
    eq2 := 2+7*x-3*y-5*x*y=0;
    solve({eq1,eq2}, {x,y});
end proc();
                                             29
                              {x = -5/3, y = --}
                                             16

Presumably in your actual problem that won't suffice, that is, you need to pass the results of one into the next.  You can do so with

proc()
local x,y,eq1,eq2,sol1,sol2;
    eq1 := 5+3*x=0;
    eq2 := 2+7*x-3*y-5*x*y=0;
    sol1 := solve(eq1, {x});
    sol2 := solve(eval(eq2, sol1), {y});
    sol1 union sol2;
end proc();
                                             29
                              {x = -5/3, y = --}
                                             16

Actually cat is the older (but better) way.  At one time the catenation operator was the dot (.).  It was changed to || to allow using the dot for Matrix/Vector operations.

No, you want the combinations, not the permutations.  Doug's algorithm will work, but instead of interpreting its output as the locations for the zeros, interpret it as the locations of the nonzeros (and adjust k accordingly). So for your example do

Q := (n,k) -> remove( L -> has( L[2..-1]-L[1..-2], 1 )
                      , combinat:-choose(n+k,k)
                     ):
Q(2,2);
               [[1, 3], [1, 4], [2, 4]]

Each list represents where the integers 1 and 2 go, the rest are zero.