Here is a much neater encoding of that theorem. In particular, pay attention to the angle-bracket matrix constructor < ...; ...>, which greatly simplifies the construction of block matrices. Also note the use of simultaneous multiple assignment (..., ...):= (..., ...), which eliminates the need for intermediary variables when updating multiple recurrences. (The parentheses are not required for multiple assignment; I just use them for clarity.)

Domineering:= proc(m::posint, n::posint, x::name, y::name)
local G0:= <<1>>, G1:= <<0>>, Z:= <<0>>;
    to n do
        (G1,G0,Z):= (<y*G0, Z; Z, Z>, <G0+G1, x*G0; G0, Z>, <Z, Z; Z, Z>)
    od;
    (G0^m)[1,1]
end proc
:                
Domineering(3,3,x,y); #Example call

The above uses 1D input (aka Maple Input).

In the 4th Reply to this Post, I wrote a tridiagonal solver that's much faster than anything offered by the LinearAlgebra package. As you suggest, my method uses prefactoring followed by back substitution. Thus, if the same coefficient matrix is used with different right-side vectors (as is typical), only the back substitution needs to be repeated.

Like this:

d:= x-> x^2+1: #denominator
Adj:= x-> 2*d(x)-x^3;
R:= 2 - (x^3+CurveFitting:-PolynomialInterpolation(
    [[1, Adj(1)], [2, Adj(2)], [3, Adj(3)]], x
))/d(x);

As is well documented, ListTools:-Classify is much more efficient than ListTools:-Categorize. It's also much more efficient than Acer's GL, and it can be used for your purpose via a two-line procedure:

GraphEigenvalueClassify:= (L::list([{Graph, name}, list(algebraic)]))->
    ([lhs,rhs]~@op@(g-> [map(`?[]`, g, [1])[]])~)(ListTools:-Classify(g-> g[2], L))
:

Here is an efficiency comparison:

Download grouplist_test.mw

You haven't indicated whether your lists of eigenvalues are always sorted. If they aren't, and you want to group by sorted eigenvalues, then simply change g[2] to sort(g[2]) in my procedure. This has neglible effect on the efficiency.

Most of the new syntax features introduced in Maple 2018 and Maple 2019, such as do-until loops, only work in 1D input (aka Maple Input) mode. If you change your input mode, then you'll have no syntax errors. The issue raised by Tom Leslie is completely different. It's a logic error---an infinite loop. They don't give error messages.

By the way, I like your neat coding style and your use of the newer syntax.

To generate a sample of size 99, for example, do

X:= Statistics:-RandomVariable(ProbabilityTable([.2, .5, .1, .1, .1])):
S:= Statistics:-Sample(X, 99):
trunc~([seq(S)])

If the support is not an initial segment of positive integers, then use EmpiricalDistribution instead of ProbabilityTable.

There's a factor of 1/sqrt(x) at the start of the expression. Thus, you need to use 

limit(asa, x= 0)

rather than

subs(x= 0, asa)

Do

remove(type, ListTools:-Flatten(w), {constant, linear})

Or, if you want to preserve the sublist structure (no matter how complicated), do

evalindets[2](w, list, type, remove, {constant, linear});

In either case, the resulting list of nonlinear terms can be counted by

nops(ListTools:-Flatten(%))

Maple has no problem doing the animation. The problem was you not stating the problem clearly. 

My animation is similar to Acer's, but I used a tangent transformation to smooth the rotation. 

restart:
Cons:= {4*x1+x2 <= 12, x1-x2 >= 2, x1 >= 0, x2 >= 0}:
Obj:= c1*x1 + 2*x2: 
LP:= proc(C1)
option remember;
    if C1::numeric then
         Optimization:-LPSolve(eval(Obj, c1=  C1), Cons, maximize)[1]
    else
        'procname'(C1)
    fi
end proc
:
plots:-animate(
    plot, 
    [
        (LP(tan(c1))-tan(c1)*x1)/2, x1= 0..4, x2= -1..2, thickness= 3, color= red,
        title= 'sprintf'("c1 = %4.2f, Z = %4.2f", tan(c1), LP(tan(c1)))
    ],
    c1= -Pi/2..Pi/2, frames= 300, paraminfo= false,
    background= plots:-inequal(Cons, x1= 0..4, x2= -1..2, nolines)
);

There's no way that you're going to get an interesting result for this, because Maple doesn't even know how to differentiate KummerU with respect to its first parameter. But if you just want a result, any result, instead of an error message, you could do

series(' 'KummerU' '(p, 3/2, t), p);

That's two set of single quotes, not one set of double quotes.

Your underlying plot command is nonsense because it specifies a plot with respect to two variables, x1 and x2. There can only be one such variable.

Your trouble has nothing to do with the background option.

You can't specify both the value of the step size h and the desired accuracy at the same time because they depend on each other.

You'll also need to specify either the number of steps or the final x-value. Usually n is the number of steps, so are you sure that you're interpretting the instructions correctly?

Your dx and point_list don't make sense because they're loop invariants: Their value doesn't change in the loop.

I do this:

solve({sec(x)=sqrt(2), x>=3*Pi/2, x<=2*Pi}, allsolutions, explicit);

This'll also work with RealDomain:-solve. IMO, it's usually better to put inequalities inside the solve rather than in an assuming clause. (And the help ?RealDomain says that assumptions don't work with that package.)

This usage of explicit is undocumented, but I've been using it since Maple 16 (and it may have existed before that). I just discovered it by experimentation. Both the keywords allsolutions and explicit are required to get the specific results that satisfy the inequalities. Taking the help page ?solve,details at face value, it makes no sense to me that these keywords would work together for a transcendental equation, but the above results show that they obviously do. 

I would avoid using RealDomain unless you're sure that you need it. Wanting only real solutions is not a good reason to use it.

You can provide fsolve with ranges within which to search, if you know some:

fsolve({r1, r2}, {xa1= 0..1, xb1= 0..1});
      {xa1 = 0.5048762298, xb1 = 0.06850883394}

This solution is returned immediately.

Any system of recurrence relations of finite order with initial conditions (regardless of linearity, homogeneity, or autonomy) is equivalent to an ODE IVP system where the first derivatives are interpretted as the discrete derivative diff(x(n), n) = x(n+1) - x(n), and this is extended to higher-order derivatives in the natural way: (D@@k)(x)(n) = (D@@(k-1))(x)(n+1) - (D@@(k-1))(x)(n). The new equations are formally called difference equations rather than differential equations. If that ODE IVP is passed to dsolve(..., numeric, ...), and solved with Euler's method with stepsize = 1, then the solutions are identical to the solutions of the recurrence relations. These solutions can be plotted with plots:-odeplot, which makes this a very convenient method.

Here's a utility to convert a system of recurrence relations to its ODE IVP equivalent:

RecToDsys:= proc(
    Sys::set(algebraic=algebraic), #recurrence relations expressed with indexed names
    V::list(name), #dependent variables
    n::name, #independent variable
    ICs::set(indexed(integer)=complexcons) #initial conditions
)
local Z,K,k,j;
    (Z,K):= (min,max)(op~(lhs~(ICs)));
    evalindets(
        eval(Sys, n= n-min(eval(op~(indets(Sys, specindex(V))), n= 0))),
        specindex(V), 
        t-> 
            local f:= op(0,t), k, K:= eval(op(t), n= 0);
            add(binomial(K,k)*diff(f(n), [n$k]), k= 0..K)
    ) union 
    eval(
       `union`(
            seq(
                {(
                    (D@@k)~(V)(Z)=~ 
                    `~`[`+`](seq((-1)^(k-j)*binomial(k,j)*~index~(V,Z+j), j= 0..k))
                )[]},
                k= 0..K
            )
        ),
        ICs
    )    
end proc
:

In order to apply this to your system, I needed to make up 6 initial conditions. You could probably come up with better initial conditions than I could. 

Sys:= {
    x[n]=1/3*(2*x[n-1]*y[n-1]+4*x[n-1]*z[n-1])+1/12*(2*x[n-2]*y[n-2]+4*x[n-2]*z[n-2]),
    y[n]=1/3*(1/4*x[n-1]*z[n-1]+y[n-1])+1/12*(1/4*x[n-2]*z[n-2]+y[n-2]),
    z[n]=1/3*(x[n-1]*z[n-1]+2*y[n-1]*z[n-1])+1/12*(x[n-2]*z[n-2]+2*y[n-2]*z[n-2])
}:
Dsys:= RecToDsys(Sys, [x,y,z], n, {x[0]=-1, x[1]=0, y[0]=1, y[1]=0, z[0]=-1, z[1]=-2});
Dsol:= dsolve(Dsys, numeric, method= classical[foreuler], stepsize= 1):
Dsol(10);
plots:-odeplot(Dsol, [x,y,z](n), n= 0..10);

For the initial conditions that I tried, the solutions either quickly blew up (diverged to infinity) or quickly converged to 0. But you may know the proper initial conditions. My procedure does NOT require that the initial conditions start at index 0.