@mthkvv As far as I know, it's not possible to programmatically set the fps (frames per second) with any of the animation commands: plots:-animate, plots:-display, or Explore. You must set it from the toolbar (or context menu) animation controls, which only appear when the animation is the current on-screen "context".

If you save the animation to a .GIF file, then you can change the fps using numerous freely downloadable GIF-editing apps.

An often-used workaround is to repeat frames. So, to change from the default 10 fps to x fps, repeat each frame round(10/x) times. However, if you do too much of this, the GUI can become overwhelmed by the amount of data and become unbearably sluggish. It largely depends on the amount of numerical data in each frame.

Change f to

f:= unapply(convert(g(x), surd), x);

I'm not 100% sure that this will work with D in Maple 13. Let me know. If it doesn't work, I can give you something else.

But regardless of whether it works with D or not, you should change f as shown. Here's a rule to help decide between x-> ... and unapply(..., x): If the right side of the arrow (in this case convert(g(x), surd)) contains a symbolic operation that doesn't depend on the specific value(s) of the variable(s) on the left side of the arrow (in this case x)---it only depends on x as a symbol---then use unapply. The convert(g(x), surd) is a purely symbolic operation: It has no dependence on the numeric value of x.

If you fail to follow this rule, the result will still often "work" in a sense, but if it does work, it'll always be extremely inefficient. The most-common violations of this rule involve putting diff or int on the right side of the arrow (in which case the result usually won't work with numeric values). (I'm not saying that you should never put those on right side of the arrow; you need to consider the situation from the symbolic perspective explained in the previous paragraph.)

This functional-compostion diagram illustrates the above. I copied it from help page ?operators,D:

       function application -->
 f   -----------------------------   f(x)
           <-- unapply
 |                                    |
 |                                  d | ^
 |                                  i | |
 |D()                               f |
 | |                                f | i
 | |                                  | n
 | V                                | | t
 |                                  V |
 |                                    |
       function application -->
D(f) ----------------------------   f'(x)
           <-- unapply

Note that you've erroneously put a space after GAMMA. It must be, for example, GAMMA(alpha+1) instead of GAMMA (alpha+1). This type of error is only possible in 2-D Input.

A vertically stacked array of three plots can be made by

restart:
F:= (c::list(numeric))-> 
    alpha-> 
        t-> add(c[i]*(t^alpha)^(i-1)/GAMMA((i-1)*alpha+1), i= 1..nops(c))
:
Alpha:= [0.2, 0.6, 0.8, 1.0]:
tr:= 0..1:
c:= Record(
    "R"= [0.46, 9.1625, 8.8318, 11.6888],
    "L"= [0.32, 0.09282, 2.1126. 3.9028],
    "T"= [0.52, 0.0569, 0.0243, 1.3102]
):
plots:-display(
    <seq(
        plot(
            F(c[x])~(Alpha), tr, legend= (alpha=~ Alpha),
            title= typeset(x, " for various ", alpha),
            labels= [t, x]
        ),
        x= exports(c)
    )>,
    view= [tr, DEFAULT]       
);

If you want all the plots to have the same scale on the vertical axis, change DEFAULT to
0..max(seq(max(F(c[x])~(Alpha)(rhs(tr))), x= exports(c)))

The plots can also be horizontally arrayed.

SumNeighborhood:= (G::GRAPHLN, v)->
    add(GraphTheory:-Degree~(G, GraphTheory:-Neighborhood(G, v)))
:

Make it a list: numelems([b]);

My procedure for this uses neither a remember table nor recursion. It does have something like a "reverse remember table". The sequence values are stored internally as hardware floats for faster computation.

This procedure is about 29 times faster than Kitonum's procedure at generating the first 2^13 entries of the sequence.

A:= proc(N::And(posint, Not({1,2})))
local 
    n, m, p, U:= table([HFloat(1)= 1]), 
    A:= rtable(1..N, [1,0], datatype= float[8])
;
    for n from 2 to N-1 do
        p:= A[n];
        m:= U[p]; 
        A[n+1]:= `if`(m::posint, n-m, min(abs~(p -~ A[..n-1])));
        U[p]:= n       
    od;
    seq(trunc(a), a= A)
end proc
:
S:= CodeTools:-Usage([A(2^13)]):
memory used=151.81MiB, alloc change=27.49MiB, 
cpu time=453.00ms, real time=460.00ms, gc time=0ns

S[..30];
[1, 0, 1, 2, 1, 2, 2, 1, 3, 1, 2, 4, 1, 3, 5, 1, 3, 3, 1, 3, 2, 
  10, 5, 8, 2, 4, 14, 4, 2, 4]

Simply add the option view= [0..24, 0..100].

This bug is very easy to work around by passing the color to display:

plots:-display(
    plots:-matrixplot(A, heights= histogram), 
    color= "X11 Thistle1"
);

There are numerous plotting commands (I've especially noticed them in package Statistics) with similar shortcomings, which may pertain to color or other options, that can be worked around in the same way. In some cases, it's necessary to also include overrideoption as an option to display.

I don't think that the Answers presented thus far can truly be called "mapping" in the fullest sense of that word because they don't produce a Matrix as output. Here's a solution that does. First, I define a procedure analogous to map:

DimMap:= proc(f, A::rtable, d::posint:= 1)
local i, D:= rtable_dims(A), R:= rtable(A);
    for i from lhs(D[d]) to rhs(D[d]) do
        R[D[..d-1], i, D[d+1..]]:= f(R[D[..d-1], i, D[d+1..]], _rest)
    od;
    R
end proc
:

The 3rd argument is the dimension over which to operate: 1 for rows, 2 for columns, and higher-dimensional Arrays are handled also.

Usage is just like map:

Sx:=1/sqrt(2)*Matrix([[0,1,0],[1,0,1],[0,1,0]]);
lam,v:=LinearAlgebra:-Eigenvectors(Sx);

DimMap(LinearAlgebra:-Normalize, v, 2, 2);

The 3rd argument, 2, is the dimension to map over. The 4th argument, 2, is the extra argument to Normalize.

Regarding your code: Matrices should be indexed with square brackets unless you're sure that you know what you're doing with the round brackets. I recall about two months ago that you had a surprising issue caused by using round brackets. I thought that you learned your lesson then. Round brackets aren't simply an alternative syntax; their semantics are very different (especially for Arrays).

Although VV's method is fast and exact, a further substantial time improvement can be made by integerizing the matrix of rational polynomials. This is done by multiplying it by its lowest common denominator. Also, the technique below can return the LU factorization of A, which'll be beneficial if you'd like to use different bs with the same A.

#This procedure returns a factorization of A, A = d*A1, where d 
#is rational and A1 is a matrix of polynomials with integer 
#coefficients.
IntPolyMatrix:= proc(A::rtable)
local
    x, 
    AR:= convert(A, rational),
    d:= ilcm(seq(x, x= denom~(AR))),
    V:= indets(A, name)
;
    1/d, rtable(d*AR, datatype= `if`(V={}, integer, polynom(integer, V)))
end proc
:
#This procedure returns a factorization of A, A = d*P.L.U, where d is
#rational and P, L, and U are the standard matrices returned by
#LUDecomposition.
LUIntPolyMatrix:= proc(A::Matrix)
local A_denom, A_int;
    (A_denom, A_int):= IntPolyMatrix(A);
    (    
        A_denom, 
        LinearAlgebra:-LUDecomposition(
            A_int, 
            ':-method'= ':-FractionFree', 
            ':-output'= [':-P', ':-L', ':-U']
        )
    )
end proc
:
#This procedure returns the exact matrix/vector X such that A.X = b.
#X is a matrix of rational functions with integer coefficients.
SolveIntPolyMatrix:= proc(A::Matrix, b::{Vector, Matrix})
local b_denom, b_int, dPLU;
    dPLU:= LUIntPolyMatrix(A);
    (b_denom, b_int):= IntPolyMatrix(b); 
    b_denom/dPLU[1]*LinearAlgebra:-LinearSolve([dPLU[2..]], b_int)
end proc
:
#Time test (using the A and b from your worksheet):
X:= CodeTools:-Usage(SolveIntPolyMatrix(A,b)):
memory used=277.44MiB, alloc change=0 bytes, 
cpu time=1.47s, real time=1.36s, gc time=250.00ms

#Verify correctness:
[seq](x, x= simplify(convert(A, rational).X - convert(b, rational)));
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 0, 0, 0, 0]

Yes, it's easy. To do it completely, there are 2-1/2 or 3 steps:

1. Use the command LibraryTools:-Save(x, y, z, ..., filename) to create a library named filename and add x, y, and z to it. The same command can be used to add new names to an existing library. Anything that can be assigned to a name can be saved; they don't need to be procedures.

2. The global variable libname determines which libraries are visible and the order in which they are accessed (very important if the same name is used in multiple libraries). libname contains the names of the directories that Maple will look in for the libraries.

3. To make it fully automatic as you describe, you must update the value of libname in your initialization file. Make sure that you update it, not just set it. For example,

libname:= libname, "the new directory":
or
libname:= "the new directory", libname:

depending on the access order you want in case of duplicate names.

I will assume that this is an initial value problem (IVP). If it's a boundary value problem (BVP), then the code below will need to be modified a bit. It's definitely an IVP or a BVP, but which one it is cannot be determined from the text snippet provided.

If you can supply numeric values for the 4 initial conditions and 2 parameters, then the rest is easy, and done below. I just made up some numbers for those 6 values. The key commands are
dsol:= dsolve(..., numeric, parameters= [...])
and
plots:-odeplot(dsol, ...)

restart
:
ODEs:= <
    #1.85:
    diff(chi(t),t$2) + 3*H*diff(chi(t),t)       
        + (exp(-sqrt(8/3)*chi(t))*phi(t) 
            - exp(-sqrt(2/3)*chi(t))*(phi(t)-zeta*diff(phi(t),t$2)^2)
          )/sqrt(24),

    #1.86:
    zeta*(diff(phi(t),t$2)+diff(phi(t),t)/3*(9*H-sqrt(6)*diff(chi(t),t)))
        - exp(-sqrt(2/3)*chi(t))/4 + 1/2
>; 
#Parameterized initial conditions:
ICs:= phi(0)=phi0, D(phi)(0)=dphi0, chi(0)=chi0, D(chi)(0)=dchi0
:
#Extra expressions:
Extra:= <
    #1.87:
    Omega__m = 
        1 -
        (2*diff(chi(t),t)^2 + 
            exp(-sqrt(2/3)*chi(t)) *
                (diff(phi(t),t)^2 + phi(t)*(1-exp(-sqrt(2/3)*chi(t))/2))
        )/12/H^2,

    #1.88:
    q = 
        -1 - 3*Omega__m/2 - diff(chi(t),t$2)^2/2
            - zeta/4*exp(-sqrt(2/3)*chi(t))*diff(phi(t),t)^2,

    #1.89:
    omega__DE = 
        (diff(chi(t),t)^2 - exp(-sqrt(2/3)*chi(t))*diff(phi(t),t)^2)
            /3/Omega__m - 1,

    #1.90:
    Omega__DE = 
        (2*diff(chi(t),t)^2 + 
            exp(-sqrt(2/3)*chi(t)) *
                (diff(phi(t),t)^2 +
                    phi(t)/H^2*(1-exp(-sqrt(2/3)*chi(t))/2)
                )
        )/12
>;      
Params:= [H, zeta, phi0, dphi0, chi0, dchi0]:
dsol:= dsolve({seq(ODEs), ICs}, numeric, parameters= Params):
#I just made up some values:
Param_vals:= Params=~ [.1, .2, .3, .4, .5, .6];
 Param_vals := [H = 0.1, zeta = 0.2, phi0 = 0.3, dphi0 = 0.4, 
   chi0 = 0.5, dchi0 = 0.6]

dsol(parameters= Param_vals):
plots:-odeplot(
    dsol, 
    eval[recurse](
        `[]`~(t, [seq](rhs~(Extra))),
        {Param_vals[], seq(Extra)}
    ), 
    t= 0..2,
    legend= [seq](lhs~(Extra))
);

It looks like you're trying to compute the stationary vector s of the Markov chain whose transition matrix is A. For a Markov chain, the matrix Id-A will always be singular. As you've noted, you need to add a side relation; in this case, the appropriate side relation is that the sum of the entries of s is 1. Once that is done, you could use a Moore-Penrose inverse to get s, but that's not the best way. 

Below, I show three ways: 1) Least squares solution of linear system, 2) Moore-Penrose, 3) Eigenvector. In this particular case of a very small A, they produce identical results (i.e., identical at the default precision Digits=10). But for larger A, I think that 2) and 3) will take more computation time and they'll have more round-off error.

In the code below, I use ^+ as the transpose operator. Depending on your Maple version and your input mode (1D or 2D), you may need to change that to ^%T.

restart;
LA:= LinearAlgebra:
# Convenient plaintext vector display format:
Ans:= (v::uneval)-> 
    printf("\t\t%a = <%9.7g>\n", v, simplify(eval(v), zero)
):

#Your transition matrix:
A:= <0.5661180126, 0.4338819876; 0.8316071431, 0.1683928571>;
n:= upperbound(A)[1]; #matrix size
#
# Method 1: Least squares solution of linear system
#
Id:= rtable(identity, (1..n)$2, subtype= Matrix);
A1:= <Id-A | <(1$n)>>^+; #include coeffs of side relation.
R:= <(0$n), 1>; #right-side vector, including side relation.
s:= LA:-LeastSquares(A1, R)^+:  Ans(s); 
             s = <0.6571429 0.3428571>

#Check:
Ans(s.A); Ans(s.A-s);
             s . A = <0.6571429 0.3428571>
             s . A-s = <1.000001e-10 1.000000e-10>

#
# Method 2: Moore-Penrose
#

#You don't need "method= pseudo" if the matrix isn't square.
s1:= (LA:-MatrixInverse(A1).R)^+:  Ans(s1);
             s1 = <0.6571429 0.3428571>

#Compare:
Ans(s1-s);
             s1-s = <1.110223e-16 -1.110223e-16>

#
# Method 3: 
#     The stationary vector is the dominant left eigenvector of the 
#     transition matrix normalized in the 1-norm
#
s2:= LA:-Normalize(LA:-Eigenvectors(A^+)[2][..,1], 1)^+:  Ans(s2);
             s2 = <0.6571429 0.3428571>

#Compare:
Ans(s2-s);
             s2-s = <-8.689216e-11 -7.542467e-12>

#That s2 contains complex values with 0 imaginary parts that could be 
#easily removed if that's desired.

If you use the options timestep and spacestep to pdsolve(..., numeric) to set those values to the same values used in your custom-coded solution, then the difference of the two curves will be much reduced. So, starting with 

pds:= pdsolve(
    PDE, IBC, numeric,
    time = t, range = 0 .. 1, timestep = 0.01, spacestep = 0.1
):

the final plot is 

You can see the contents as a list by

[seq](u[1..11, 1]);

You can plot it by

plot(<<seq(1..11)> | u[1..11, 1]>);