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;

As Kitonum points out, this is done by the command GraphTheory:-SequenceGraph. However, this will return a graph with the vertices generically labeled from 1 to n. I got the impression from your Question that you wanted to use your own vertex labels. So here is a procedure to do that. The procedure takes two arguments. The first is the list of degrees. The second argument is optional and is the corresponding list of vertices. If omitted, generic 1-to-n labels will be used. The procedure returns FAIL if the degree sequence can't possibly be made into a graph.

RealizeDegreeSequence:= proc(
     Degrees::list(nonnegint),
     Vertices::list({symbol,string,integer}):= [$1..nops(Degrees)]
)
local
     Edges:= table(), e:= 0,
     D:= Vector(Degrees),
     V:= Vector(Vertices),
     P, k
;
     if max(D) >= op(1,D) then  return FAIL  end if;
     do
          (D,P):= sort(D, 'output'= ['sorted', 'permutation']);
          V:= V[P];
          if D[-1] = 0 then  break  end if;
          for k from 2 to D[-1]+1 do
               if D[-k] = 0 then  return FAIL  end if;
               D[-k]:= D[-k] - 1;
               e:= e+1;
               Edges[e]:= {V[-1],V[-k]}
          end do;
          D[-1]:= 0
     end do;
     GraphTheory:-Graph(Vertices, convert(Edges, set))
end proc:

I'd rather do the example with Sum rather than sum because the evaluation of the sums gives results that obsfucate the example, but I wrote the code so that it works for either. Using subsindets, it's a short one-liner.

Download subsindets.mw

I found two substantial differences between your Mathematica code and your Maple code (from your PDF versions). In the Mathematica code, the initial condition for h is

G:= 0.005*exp(-(x-0.7)^2/(1/500)) + hin;

In the Maple code it's

F:= 0.01*exp(-(x-0.4)^2/(1/300)) + hin;

The other difference is that in the Mathematica code, the upper bound of x is w = 1, whereas in the Maple it's l = 2.

If I use your Mathematica conditions in the Maple code, the animation runs much faster.

I redid your worksheet in a much cleaner style (even though it was already quite well organized), which is much less prone to user errors. All floating-point constants are in a list of parameters that is not applied to the system until needed. Thus you can see the PDEs in their fully symbolic form. I changed zita to zeta[0] and gama to gamma (and made gamma local) so that these would print nicer. 

Then I used numeric differentiation (fdiff) applied to the pdsolve solution to compute the residuals of the two PDEs over a grid of 18 interior points. The residuals seem reasonably small to me, except near (x,t) = (1,5). Also, read the help page ?pdsolve,numeric,errorcontrol.

Download PDE_residuals.mw

The Maple command for (infinite) series is sum:

sum(f, n= 2..infinity);

Maple returns infinity, meaning that the series diverges. To restrict to odd terms, do

sum(eval(f, n= 2*k+1), k= 1..infinity);

If you actually mean the convergence or divergence of the sequence rather than the series, use the limit command:

limit(f, n= infinity);

which returns 0.

Considering the triviality of this problem, I suggest that you look up the Limit Comparison Test in your textbook or on the Web. If you knew this test, you would be able to look at this problem and instantly know that the series diverges.

This Answer is to answer the OP's followup question requesting a procedure that determines whether an arbitrary polynomial can be produced as a linear combination of a set of polynomials, and, if so, returning the coefficients.

I decided that it would be better to just write a separate procedure for this rather than rolling it into my procedure LinearCombo from an earlier thread. So I wrote a little procedure PolyLinearCombo that calls my version of Kitonum's CoeffsOfLinComb. I changed the declaration-of-variables argument to be optional; by default, all nonconstant names will be considered variables. If you want to use parameters, then include the variables argument as a set.

For various reasons, it's better that a procedure always return the same number of return values. So, this procedure returns true and a list of coefficients if the second argument can be produced as a linear combination of the elements of the first argument, and it returns false, [] if there's no such combination.

I decided to disallow 0 as the second argument. If you want nontrivial coefficients to produce 0, then just call CoeffsOfLinComb.

CoeffsOfLinComb:= proc(
     L::list(polynom),
     V::set(name):= indets(L, And(name, Not(constant)))
)
local
     c, k, C:= {c[k] $ k= 1..nops(L)},
     S:= solve({coeffs(expand(`+`((C*~L)[])), V)}, C),
     F:= indets(rhs~(S), name) intersect C=~ 1
;
     eval([C[]], eval(S,F) union F)
end proc:

PolyLinearCombo:= proc(
     F::list(polynom),
     f::And(polynom, Not(identical(0))),
     V::set(name):= indets([F,f], And(name, Not(constant)))
)
local C:= CoeffsOfLinComb([f, F[]], V);
     if C[1]=0 then (false, []) else (true, -C[2..] /~ C[1]) end if
end proc:

Examples of use:

Dependent example:

F:= ['randpoly([x,y])' $ 22]:
f:= randpoly([x,y]):
C:= PolyLinearCombo(F, f, {x,y});

Independent example:

F:= ['randpoly([x,y])' $ 22]:
f:= randpoly([z,y]):
C:= PolyLinearCombo(F, f, {x,y,z});

The above example becomes dependent if z is treated as a parameter:

F:= ['randpoly([x,y])' $ 22]:
f:= randpoly([z,y]):
C:= PolyLinearCombo(F, f, {x,y});

I'd condsider this a bug: The value of Digits and other environment variables (which have been modified from their defaults) should be transcribed to all nodes. Another workaround is to put

Digits:= 30;

into your initialization file.

I agree with Roman that you should run the evalf also in parallel; timewise, it is a significant part of the computation. But I do not agree that you should necessarily use numeric integration. So, change the Mapped procedure from int to evalf@Int or evalf@int depending on whether you want to use numeric or symbolic integration.

Your problem is that pdf1(i), etc., will evaluate to a numeric value even if i is symbolic (i.e., even if i has no numeric value). There are two ways to correct this. The preferred way is to make the following the first executable line of pdf1:

if not x::numeric then return 'procname'(args) end if;

The crude way is to put unevaluation quotes around the call to pdf1:

Sum('pdf1'(i), i= 3..3),

etc.

There is no bug.