It is not possible to use plots:-animate because HeatMap is implemented to display a background picture.
Explore can be used, but I am not sure whether it works in Maple 2016. In Maple 2017+ it's ok.

restart;
p := k -> Matrix(3, 3, (i,j) -> (i+j+k) mod 3):
f:=proc(k) Threads:-Sleep(0.2): Statistics:-HeatMap(p(k)) end:
Explore(  f(k), k=0..8, animate,loop, autorun);

When a variable is implicitely declared local, a warning message appears.
The rules for this are simple, see ?local

In

foo:= proc() 
  local x;
  plot(sin(x),x=-Pi..Pi); 
end proc:

local is indeed needed, otherwise x:=10  at top level will produce an error.

Alternatively you may use
plot(sin('x'),'x'=-Pi..Pi);

All the floating point computations are done with 1 significant digit. So, sqrt(5) is approximated to 2.  and  (2. -1.)/2.  is 0.5.

See also:

restart;
Digits:=1;
                          Digits := 1
a:=[seq(i/10., i=1..9)];
       a := [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
b:=[seq(i/10., i=9..1,-1)];
       b := [0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]
add(a) = add(b);
                            6. = 4.
evalf[2](add(a)) = evalf[2](add(b));
                           4.5 = 4.5

Maple (at least the recent versions) can compute directly the double integral (without passing to polar coordinates):

int(sqrt(x^2+y^2), x = 0 .. B/2, y = 0 .. b/2) assuming b>0, B>0;

S:=simplify(eval(SA/t, [x=y*t, Log[t]=1/L])):
collect(S,y, expand);

Now you have a polynomial in the variables (L,y) with rational coefficients and you want to guess its coefficients or a generating function. Why do you think that a simple generating function does exist? A more realistic approach would be to come back to the original problem (which produced the expression) and try there.

.

n >= ceil(fsolve(1/(n+1)! = 0.00001));
                             8 <= n
n >= ceil(fsolve(1/(n+1)! * exp(0.1) = 0.00001));
                             8 <= n

Seems to be a bug. Workaround:

ex:=D[2](eta)(t,x)+D[2](phi)(t,x,0)+D[1](phi)(t,x,0);

inds:=[indets(ex)[]]:  cinds:=map(convert, inds, Diff):
subs(inds=~cinds, ex);

Download Interpolate.mw

a:=sqrt(6)/6 + sin(5/6) + exp(6/7) + 567:

f := u ->  `/`(subs(6=b, [op(u)])[]):
subsindets(subs(6=b,a), fraction, f);

For non-numeric coefficients you may use
signum(lcoeff(expr));

@AliahNiu 

For such functions Maple cannot help much, we must use some maths.
You have two (possible) singularities at x=0 and x=Pi.
For the integral to exist near 0 a necessary and sufficient condition is  (a+1/2)*p+1-d < 0.
Near Pi the condition is (b+1/2)*p+1-d < 0.
So, the integral exists iff both conditions hold.

Note. Here the integral was considered a an improper Riemann one (in this case this being equivalent with Lebesgue integrability). If you want standard Riemann integrability (equivalent here with continuity) the conditions are similar but more restrictive, namely <= -1  instead of  < 0.

Download check-pde.mw

Edit. Actually for u() given by Maple it is easy to compute  u(1, 1/2)  = 1/12, so definitively wrong!

I'll choose the space for you.

u1 := (x, y)-> 1-exp(a*x)*cos(2*Pi*y):
u2 := (x, y)-> a*exp(a*x)*cos(2*Pi*y):
Norm2 := u -> int(u^2, 0..1,0..1)^(1/2):
simplify( Norm2(u2)/Norm2(u1) );

Create in a folder e.g.  d:/temp   a file named eps.mpl   containing:

with(plots):
plotsetup(cps,plotoutput="d:/temp/ex.eps",plotoptions="portrait,noborder");
densityplot(x,x=0..1,y=-0..0.1,colorstyle=HUE,style=patchnogrid,scaling=constrained,view=[0..1,0..0.1],axes=none);
# plotsetup(default);

Now start the Command-line Maple and type

read "d:/temp/eps.mpl";

You will find the desired file "d:/temp/ex.eps"

P.S. This workaround uses the fact that the graphics 2D eps drivers for the "classic" Maple work, unlike the "standard" ones.
Unfortunately, for 3D, this does not work, the "classic" 3d drivers are wrong (at least in my computer).
I cannot understand why Maplesoft could not find a solution for this very old problem. When I need a 3d Vector graphics file, I have to  use Maxima because Maple cannot obtain a correct one.

I wrote this version a few days ago. It's af course similar to Carl's evalhf version, but >7 times faster due to compilation.
It is a pity that the compiler does not accept sub-procedures; we would have then a really useful tool to create fast numeric programs.

Sj:=proc(d::integer[8], j::integer[8])::float[8];
local k::integer[8], r::float[8],  
      x::integer[8],y::integer[8],n::integer[8], P::integer[8], u::integer[8], v::float[8];
option autocompile;
r:=0.0;
for k from 0 to d do
   n:=d-k; x:=16; y:=1; P:=8*k+j;
   if n = 0 then u:=1
   else
      while n > 1 do
         if type(n,even) then x := x^2 mod P; n := n/2;
         else  y := x*y mod P; x := x^2 mod P; n := (n-1)/2 fi;
      od;
      u:=x * y mod P
   fi;
   r:=r+u/evalf(8*k+j);
od;
r:=frac(r);
v:=1./16.;
for k to 14 do
  r:=r+v/evalhf(8*(k+d)+j);
  v:=v/16;
od;
frac(r);
end:

BBP:= proc(d::posint)
#d is the starting digit position. 9 hex digits will be printed
option hfloat;
local dm1:= d-1;
Digits:= trunc(evalhf(Digits)) + 3;
round(16^9*frac(4*Sj(dm1,1) - 2*Sj(dm1,4) - Sj(dm1,5) - Sj(dm1,6)+5));
nprintf("%09X",%); #convert(%,hex);
end proc:

BBP(1000000);
     26C65EC54