These numbers can be extracted very simply and quickly as map2(-op, 1, roots(numer(q))) and map2(-op, 1, roots(denom(q))) or a slight variation of those if you need the multiplicities and/or assignment to a[1], a[2], ..., b[1], b[2], .... 

To plot column Vector A as a function of same-length column Vector X, simply do

plot(< X | A >);

I don't mean to disparage Mac Dude's Answer; certainly it is what one would use in most cases. However, it is so obvious that I'd be surprised if it hadn't already occurred to the OP. Thus, I wonder whether the OP had something else in mind: to express the iteration over two variables with a (single) loop. This is possible, and it can be done for an arbitrary number of variables.

The easiest case is when the set of desired n-tuples can be expressed as a Cartesian product A^n for some set A. In the example that the OP gives, this would be {$1..10}^2. For these cases, combinat:-permute can be used. Like this:

for ij in combinat:-permute(op~([[$1..10] $ 2]), 2) do
     (i,j):= ij[];
     #rest of loop
end do:

Cases where the n-tuples are the Cartesian product of different sets can also be handled. This requires an auxillary procedure---a substitute for combinat:-permute---that returns the Cartesian product. Here is one such:

CartProdSeq := proc(L::seq(list))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
eval([subs(Seq=seq, foldl(Seq
                             , [cat(i,1..nargs)]
                             , seq(cat(i,j)=L[j],j=nargs..1,-1)
                            ))]);
end proc:

And an example of using that in a loop: Suppose we wish i to vary from 1 to 10, j from 1 to 9, and k from 1 to 8:

for ijk in CartProdSeq([$1..10], [$1..9], [$1..8]) do
     (i,j,k):= ijk[];
       #rest of loop
end do:

MDD, 

Please test the procedures in the attached worksheet. The first procedure, Homogenize, does what you want to the list of polynomials, F. The second procedure, MDD, is just a trivial one-liner that allows for your distinguished single polynomial, f.

Download Homogenize.mw

Symbolically eliminating the parameter is only possible in very special cases. (Whether this is such a case I don't wish to discuss in this Answer.) Cases where Y can be completely expressed as a function of X are so trivial that they rarely are used in practice. One obstacle is that a range of the parameter must be specified so that Y is indeed a function of X rather than there being possibly multiple Y for each X. However, if you do specify such a range of the parameter, then it may be possible to at least give a numerical procedure that expresses Y as a function of X. Here is one for your example:

X:= (theta)-> cos(theta) + 0.8e-1*cos(3.*theta):
Y:= (theta)-> -sin(theta)+ 0.8e-1*sin(3.*theta):

YY:= proc(XX::numeric, theta_range::range)
local theta;
     Y(fsolve(X(theta)=XX, theta= theta_range))
end proc:

plot('YY(XX, 0..Pi)', XX= -1..1);

This is crude, and can be improved with some error checking, such as making sure that the fsolve actually returns a numeric value. 

To correct the "too many levels of recursion" error, in the definition of function g0, change sum to add. The problem with sum is that it tries to simplify these horrendously long nested piecewise expressions using symbolic j. If you use add, then it only uses definite integer values of j, and the piecewise expressions can be immediately resolved.

Apply CylindricalAlgebraicDecompose separately to the numerator and denominator of your rational function. You can use numer(r) and denom(r) to extract the numerator and denominator of rational function r. In some cases, it may be beneficial to first use normal(...), as in denom(normal(r)). See ?numer for some examples.

If you don't want to go the simplify(..., symbolic) route---which can be risky---it's easy to have a one-liner procedure that extracts all the bases of `^` expressions and and sets them  > 0. 

Download AllBasesPositive.mw

Here's how to fix this error. I'll assume that you still have Maple 16 installed. If you don't, let me know, and I'll post here the necessary code.

There are three procedures necessary for your program that exist in Maple 16 that are missing in later Maple. Their names are zero_one, polyvariations, and midpoint. These procedures are defined when you execute the (deprecated) command readlib(realroot). We need to extract the code of these procedures from Maple 16 and put it in your code file. Here's how.

1. In a Maple 16 worksheet, issue the command

readlib(realroot):

2. Issue the command

showstat~([zero_one, polyvariations, midpoint]);

3. Using your favorite text editor (I used Microsoft Windows Notepad), cut-and-paste the code of the three procedures to the end of your code file proving.txt.

4. Using the text editor, remove the line numbers from the procedures, being very careful that you don't remove anything else. It's okay to remove or add spaces at the beginnings of lines.

5. Add a comment to your code about where these three procedures came from. Mention how they were extracted with readlib.

6. Save the file proving.txt.

7. Run your program in Maple 2015.

I have done this, and it works.

There may still be use instances of your program that lead to similar errors, or will eventually lead to similar errors in future editions of Maple. To help in finding those errors now rather than when those future editions come, I recommend that you comment out all uses of the command readlib and then do some rigorous testing of your program---something that uses all of its procedures.

Several changes were necessary to get this to work, all of which I point out in my modified worksheet below. Most significant is turning off hardware floats and specifying the inflection points of the distribution (see ?Statistics,Sample). It's a shame---an unnecessary limitation of Statistics:-Sample---that it does not automatically revert to software floating point if the hardware floating point doesn't work.

I redid your worksheet using 1D input because I can't tolerate the 2D input. I put my name in front of all my comments.

Download RandomVariableSample(2).mw

How about

M:= k-> Matrix(10, 10, [1$k, 0$10-k], scan= diagonal):
M(5);
'M(k)' $ 'k'= 1..10;

The evalf[n](expr) doesn't mean "give me the value of expr accurate to n digits". Many users have that wrong impression. It means (approximately) that at least n digits are used in all the subcomputations. That "give me the value..." thing is practically impossible to implement. For the reasons, see the first chapter of almost any textbook on numeric computation.

The code that you're running here is obviously intended to break or to highlight the deficiencies of floating-point computation systems.

As to your second question: Many Maple commands "remember" the values from earlier computations. So, if it knows that it already computed the desired result to a higher precision, it just re-uses that. The commands restart and forget make it forget what it has remembered (I recommend the restart).

Three points:

1. Don't set Digits as low as 5 unless that's necessary to get an answer. If you only want 5 digits, just ignore the other digits. (There's also a command interface(displayprecision). The results of using that are visually hideous.) Stick with the default Digits = 10, or use my favorite "sweet spot", Digits = 15.

2. The command pade is in the package numapprox. To access it you need to refer to it as numapprox:-pade or numapprox[pade] or first issue the command with(numapprox). I strongly prefer the first of these three options.

3. solve will spend a long time trying to find exact solutions to this pair of high-degree polynomials---a likely impossible task. Use fsolve instead of solve to get decimal solutions. This only takes a few seconds. This will only give one solution. If you need more than one, let me know: For polynomial systems, it's usually easy to get all the solutions within any real bounding box.

So, change the last part of your your code to

fsolve({
     limit(numapprox:-pade(diff(f,x), x, [4,4]), x= infinity) = 1,
     limit(numapprox:-pade(t, x, [4,4]), x= infinity) = 0.
     }, {A,B}
);

     {A = -0.680502612209264e-1, B = -3.43533296567616}

Here's an even shorter and faster variant for returning the coefficient of particular term.

P:=a*x^2*y*z+b*x*y*z+c*x^2*y*z+d*x*y^2*z+e*x*y*z:
L:=[coeffs(P, [x,y,z], 't')]:

H:= table([t]=~L):
H[x^2*y*z];

       a + c

How about (assuming that alpha has a numeric value)

fracdiff(t^beta, t, alpha);

or, if alpha is symbolic, something like

fracdiff(t^beta, t, alpha) assuming alpha > 0, alpha < 1;