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

An interesting variation is to consider the "natural" continuous extension of your function obtained by substituting cos(n*Pi) for (-1)^n in the closed-form discrete solution:

combine(rsolve({B(n)=(B(n-1)+B(n-2))/2, B(0)=0, B(1)=1}, B(n)));
b:= unapply(subs((-1)^n= cos(n*Pi), %), n);
plot(b, 0..10);

@vv Okay, replace simplify(..., symbolic) with op@@2:

f:= z-> ln(exp(1/z - 1/sqrt(z^2))):  r:= t-> exp(I*t):
IntegrationTools:-Change(Int(f(z), z= r(-%Pi)..r(%Pi)), z= r(t));
value(subsindets(%, ln(exp(anything)), op@@2));