It works here.  A better technique is to call solve with the variable enclosed in a set, that returns an expression sequence of sets of equations defining the solution (rather than an expression sequence of values).  That is, do

k1 := -(ln(q/(-1+q+exp(phi)))+phi)/phi:
eqn1 := (1-q)/q = mu:
sol := solve(eqn1,{q});
                                            1
                              sol := {q = ------}
                                          1 + mu

subs(sol,k1);
                                      1
                   ln(---------------------------------) + phi
                               /       1              \
                      (1 + mu) |-1 + ------ + exp(phi)|
                               \     1 + mu           /
                 - -------------------------------------------
                                       phi

You can use LinearAlgebra:-RowDimension (and ColumnDimension). Or LinearAlgebra:-Dimension for both.

Use simplify with side-relations.  Here is a simple procedure, Identify, to create the identifications:

Identify := proc( s :: set )
local i,j,n;
    n := nops(s);
    {seq(seq(s[i]=s[j], j=i+1..n), i=1..n-1)}
end proc:

p := (x+y+z)^3:
simplify(p, Identify({x^2*y, x^2*z, y^2*x}));

             2          2      2          2              3    3    3
          9 x  z + 3 x z  + 3 y  z + 3 y z  + 6 x y z + x  + y  + z

You could specify an additional variable to algsubs, that works if you are replacing n, however, you need to remove the sqrt,

algsubs((Variance(X3))=sigma^2,k3,n);

The Maple ^ operator is syntactically not associative (more precisely, non-associative), so x^y^z is ambiguous and raises an error.  In 2D math things are different in that the action of typing the caret (^) switches to superscript mode. The switch to superscript make it natural to make the operation right-associative.

Interpreting x^y^z as (x^y)^z is non-standard.  In 2D mode one can already get that effect by typing, in 2D, x^y*z, which creates x^(y*z), which is equivalent to (x^y)^z. 

I suspect that your interpretation comes from calculator usage, where typing x^y^z computes (x^y)^z.

Use printf:

printf("Stationary point %a is a %a\n", s[i], sptype);

You might consider using the inert Power function.  It is already known to Maple (used with mod).  Here is an extension to combine that does what you want

`combine/Power` := proc(expr1)
local pwrs, rest, p, L, b, e, expr;
    if expr1 :: `*` then
        expr := map(procname,expr1);
        (pwrs,rest) := selectremove(type,expr,'Power(anything,anything)');
        pwrs := `if`(pwrs :: `*`
                     , [op(pwrs)]
                     , [pwrs]
                    );
        pwrs := [ListTools:-Categorize((a,b) -> op(1,a)=op(1,b), pwrs)];
        return rest*mul('Power'(op([1,1],L), add(op(2,p), p in L)), L in pwrs);
    elif expr1 :: `+` then
        return map(procname, expr1);
    elif expr1 :: 'Power(anything,anything)^anything' then
        expr := map(procname, expr1);
        (b,e) := op(expr);
        return Power(op(1,b), op(2,b)*e);
    else
        return expr1;
    end if;
end proc:


combine(23*Power(a,b)^2 / Power(a,c) * Power(b,c));
                                       23*Power(a,2*b-c)*Power(b,c)

Be aware that the ?GraphTheory package does not permit graphs with self-loops.  A self-loop is implied if there are any non-zero elements on the diagonal of the adjacency graph.  For example:

GraphTheory:-Graph(Matrix(2, [[1,1],[1,0]]));
Error, (in GraphTheory:-Graph) an undirected graph cannot contain loops (loop
detected at vertex 1)

What version of Maple are you using, and is it on a 32 or 64 bit machine?  This works for me in in cmaple, the Standard GUI with both 1D and 2D input.

Your description is not precise.  The term 'function' in Maple specifically refers to an unevaluated function call, say f(x[1],3).   Is that what you mean?  If so, then where are the expressions of x[1] for which you want coefficients?  Inside the a, b, and c?   While a general description of the problem doesn't hurt, a specific example is usually easier to understand and allows one to test.
 

Maple tables are hash tables. 

A straightforward way is to use a table to store the data, then convert the table to the form you need.  For example:

temp := proc(xf, dx, dt)
local cnt, Data, x, t;
    Data := table();
    t := 0;
    x := 0;
    cnt := 1;
    Data[cnt] := [t,x];
    while x < xf do
        x := x+dx;
        t := t+dt;
        cnt := cnt+1;
        Data[cnt] := [t,x];
    end do;
    return [seq(Data[cnt], cnt=1..cnt)];  # the use of cnt as the index is a trifle tricky
end proc:

Data := temp(3, 0.1, 0.2);

plot(Data);

You can change the error tolerance using the abserr and relerr parameters.  You could also use a different integration method.

This is the chain-rule, using Leibnitz notation.  Because there is some ambiguity in the expression, Maple uses a different notation:

(**) z := x(y(t));
                                                              z := x(y(t))

(**) diff(z,t);
                                                                     /d      \
                                                          D(x)(y(t)) |-- y(t)|
                                                                     \dt     /

# Here we convert D(x)(y(t)) to the diff notation,
(**) map(convert,%,diff);
                                                    / d       \|          /d      \
                                                    |--- x(t1)||          |-- y(t)|
                                                    \dt1      /|t1 = y(t) \dt     /

You can use  multiline comments:

(*
   These lines won't be executed
*)

Alas, they are not available in 2D mode.