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));

You can upload a worksheet. The error pertains to the inline display of the worksheet, not to its uploading.

Both the Answers by Rouben and by Tom assume that you're trying to solve equation p for u and then verify that solution by substitution and simplification to 0. However, my reading of your worksheet suggests that you want to substitute an expression for u into p and solve the resulting equation for q (and verify). So, which is it?

My code is for your example only. Like Kitonum, I don't understand what you want for the general case. Also, it's not clear whether you want to construct the sets (list their elements) or define them abstractly (as in a membership test).

@vv The simplify(..., symbolic) is only applied to the limits of integration (see the second argument of the subsindets), so it doesn't touch sqrt(z^2). Using that f, the integral returns 2*I*Pi (after using evalf).

@Kitonum An easy generalization to any number of dimensions, with the input and output both being vectors:

PV:= proc(v::Vector) 
local x:= `tools/genglobal`(:-_x), pv:= Vector(upperbound(v), symbol= x);
    eval(pv, solve(v.pv, {seq(pv)}))
end proc: 

Or, a simpler version, allowing the user to specify the free global stem:

PV:= proc(v::Vector, free::name= :-_X) 
local pv:= Vector(upperbound(v), symbol= free);
    eval(pv, solve(v.pv, {seq(pv)}))
end proc: 

@vv Thanks for finding that. Here's a correction:

f:= z-> z/abs(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)), simplify, symbolic));

The above returns 0. Changing f to z-> 1/z, it returns 2*I*Pi.

@Christopher2222 I agree with your policy regarding voting. For almost any legible and legitimate math/Maple Question, if the author has reputation < 11, I give it a vote up. The 11 is the cut-off below which the Question gets the red spam flag.