I said that it would be possible to construct the polynomial (G0^m)[1,1] by modular imaging and polynomial interpolation. Here I present code that does this. The problem is that this code is slower than direct computation by powering the matrix of polynomials. But, I present it in the hope that someone may be able to suggest a way to make it faster.

Domineering:= proc(
    m::posint, n::posint, x::algebraic, y::algebraic, 
    {method::identical(naive, modular):= naive}
)
uses LAM= LinearAlgebra:-Modular, CF= CurveFitting;
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;

    # Naive  method
    if [x,y]::[name$2] implies method=naive then return (G0^m)[1,1] fi;

    #-----------------------
    # Modular Imaging Method
    #-----------------------

    local 
        #-d..d will be the interpolation points (based on degree of result in x and y).
        d:= ceil(m*iquo(n,2)/2),

        #Need primes whose product is greater than the maximum possible polynomial
        #evaluation in the [1,1] position of the mth power.
        p:= 2^25, #upperbound of float[8] primes
        bits:= ilog2(p)-1, `2^n`:= 2^n, `2^m/2`:= 2^m/2,
        primes:= [
            to ceil(ilog2((`2^n`*eval(G0[1,1], [x,y]=~ d))^`2^m/2`/`2^n`)/bits) do 
                p:= prevprime(p) 
            od
        ],
        primes1, i, j, G0ij, Pxy:= table(), #integer interpolation points
        st:= time(), ps:= 0 #efficiency info (userinfo only)
    ;
    :-`mod`:= mods; #for chrem
    userinfo(
        1, Domineering, 
        "Using d =", d, ",", (2*d+1)^2, "evaluation points, and",
        nops(primes), "primes max." 
    ); 
    #
    #This loop is what I want to make more efficient. Everything else is fine.
    #
    for i from -d to d do  for j from -d to d do
        G0ij:= eval(G0, [x,y]=~ [i,j]);
        primes1:= primes[
            ..ceil(ilog2((`2^n`*max(2,abs(G0ij[1,1])))^`2^m/2`/`2^n`)/bits)
        ];         
        userinfo(1, Domineering, NoName, proc() ps+= nops(primes1); [][] end proc());
        userinfo(2, Domineering, 'i'=i, j='j', 'nops(primes1)'=nops(primes1));
        Pxy[i,j]:= chrem(
            [seq](
                trunc(LAM:-MatrixPower(p, LAM:-Mod(p, G0ij, float[8]), m)[1,1]), 
                p= primes1 
            ),
            primes1
        )
    od od;
    userinfo(
        1, Domineering, 
        "Main loop cputime =", time()-st, "with", ps, "evaluations."
    );
    
    local Px:= table(); #interpolants wrt y
    for i from -d to d do
        Px[i]:= CF:-PolynomialInterpolation([$-d..d], [seq(Pxy[i,j], j= -d..d)], y)
    od;

    #final interpolation, wrt x
    CF:-PolynomialInterpolation([$-d..d], [seq(Px[i], i= -d..d)], x)
end proc
:         
#Example usage:
p1:= CodeTools:-Usage(Domineering(6,6,x,y));
memory used=51.45MiB, alloc change=0 bytes, cpu time=235.00ms, 
real time=237.00ms, gc time=0ns
                 [long polynomial omitted]
infolevel[Domineering]:= 1:
p2:= Domineering(6,6,x,y, method= modular);
Domineering: Using d = 9 , 361 evaluation points, and 22 primes max.
Domineering: Main loop cputime = 22.078 with 5757 evaluations.
                 [long polynomial omitted]
simplify(p1 - p2); 
                 0

@acer At first, I used the approach that you suggest. The reason that I went with op(0,a) is that it leaves only one occurence of A in the whole code, making it easier to change that to another name. Of course, the name could be made a procedure parameter also.

Yes, I realize that your way is better if we're just dealing with a list of As. I mentioned that a few Replies back in the last paragraph. 

@emendes 

The type specindex(A) means any indexed name that starts with A. The spec is short for specific. 

Let a be an indexed name, for example, A[ j1, j2 ]. Then op(0,a) is A, op(1,a) is j1, and op(2,a) is j2.

@emendes I was hoping for the desired Maple output, but it's a moot point now.

@emendes Here's a simplifcation. This works because the operation on the 1st index commutes with the operation on the 2nd index.

(T1,T2):= (table([2,3]=~ [3,2]), table([2,3,5,6,7,9]=~ [3,2,6,5,9,7])):
T:= (T::table, x)-> `if`(assigned(T[x]), T[x], x):
subsindets(varA, specindex(A), a-> op(0,a)[T(T1, op(1,a)), T(T2, op(2,a))]);

This is also easier to understand, I think. The 2nd line says that the transformation is the identity except for indices listed in the respective table.

If your actual-use cases are simply lists of names that all begin with A, then you might as well use Acer's solution. I was assuming that the actual-use cases would have the A names mixed into other expressions. In that case, you need subsindets or evalindets.

@santamartina If mapping over a list, there will be one result for every member of the list.

Your code could be replaced by this recursive procedure:

p:= (A::Array, x::posint)-> 
    ArrayTools:-Append(`if`(x > 1, thisproc(A, x-1), A), A[-1]+1)
:
p(Array([4]), 3);

@emendes There's no need for an intermediary transfer to B. The whole job can be done by

(T1,T2):= (table([2= 3, 3= 2]), table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7])):
subsindets(
    subsindets(
        varA, 
        `?[]`~(A, `[]`~(anything, {indices(T2, 'nolist')})), 
         a-> A[op(1,a), T2[op(2,a)]]
    ), 
    `?[]`~(A, `[]`~({indices(T1, 'nolist')}, anything)),
    a-> A[T1[op(1,a)], op(2,a)]
);

You'll need to show me what the desired output is.

@darcd An assuming clause is only in effect for the command that it's attached to. The simplify is irrelevant. Replace simplify with combine.

@Joachim Sand Here's a pedagogical technique with which you might make progress on this: Suppose that you were designing a synthesizer for just two frequencies instead of seven. Give us the equations for that, and let's see if they're solvable. If they are, then we'll go to three frequencies. If they aren't, it'll be easier to find the error in the smaller system of equations.

@Joachim Sand Looks like you're designing a synthesizer.

I think that your units don't balance: Assuming that all of the Rs ({R__4, R__2, R__3, R__C, R}) are resistances, C is capacitance, V__i is voltage, 440 is frequency, and 1, 3, and all the powers of 2 are dimensionless; then the left-side dimension of each equation is voltage, but the right-sides are dimensionless.

@acer Great idea with the undefined. Vote up. In the case of using plot, Float(undefined) can be abbreviated to undefined.

Acer is right; I stand corrected. Vote up.

I wrote the above code for Maple 18, as claimed by the OP's header. In more recent versions of Maple, it can be replaced by the slightly more readable 

(add/numelems)~(
    map2(
        index, d, 
        remove(`=`, [ListTools:-Split(j-> data[j]>0.01, [$1..numelems(data)])], [])
    )
);

Your data only have two points in time---2013 and 2014. That's hardly enough for a time series analysis.

Also, only one number is given for the "total population" for each region. Thus, we must assume that those totals remain the same from 2013 to 2014. That hardly seems realistic.

The only thing that might be possible to assess statistically is whether the 2014 numbers are significantly different from the 2013 numbers. Even that's rather shaky because of changing population, as noted above.