@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.

@awass The standard order of evaluation is as you describe: from the inside to the outside. The vast vast majority of Maple commands evaluate this way. A few commands necessarily must operate differently, and evalf is one of those.

You can nest evalfs with different precisions. Note the difference in the following two evalf commands:

a:= 1:
evalf[5]((a+1e-5)^100);
             1.0000

evalf[5](evalf[10]((a+1e-5)^100));
             1.0010

@mmcdara The with command should never be used inside a procedure. In addition to being very time consuming, it's very unreliable. This is described on its help page. Instead, use uses LinearAlgebra.

Please stop posting things that should be considered Answers as Replies.

@RohanKarthik 

Note that I made a slight change in the code in my Answer. I think that the new version is slightly more efficient. In any case, it's easier to understand.

Some observations on the sparsity, without any conclusions yet:

1. The density (the fraction of nonzero elements) of G_(0,n) is exactly (3/4)^n. Hence, density --> 0 as n --> infinity. G_(1,n) has roughly half the density of G_(0,n).
2. The sparsity pattern of G_(0,n) is a bit unusual: It's upper triangular with respect to the non-principal diagonal.
3. Your observation is correct, but your suggestion is ineffective: Matrix A^m is not computed as A^(m-1).A; rather repeated squaring is used. For example, to compute A^9: A^2 <-- A.A, A^4 <-- A^2.A^2, A^5 <-- A^4.A, A^9 <-- A^5.A^4. Thus, only 4 matrix multiplications are needed to compute A^9.
4. The entries of G_(0,n) (not G_(0,n)^m) are closely related to Fibonacci numbers.
5. It's possible to construct polynomials by performing computations over prime-number fields and interpolating. Maple's dense matrix computations over prime fields is extremely fast.

@vv The point is to start with a form of the problem as close as possible to the original and have Maple do the rest. I think that the form that I have is as close as one can get to the original without a lengthy procedure.