A faster version.

F := proc(N::posint)
  local V:=Vector(4), p:=2;
  while (p:=nextprime(p))<N do V[(irem(p,8)+1)/2]++ od;
  [seq(V)]
end proc:

A := Matrix(2, 3, [[1, 2, 3], [3, 1, 2]]); B := Matrix(2, 3, [[1, x, 3], [y, 1, z]]);
solve(Equate(A,B), {x,y,z});

Note that solve needs a list or set for its first argument, so solve(A=B, {x,y,z}); will not work.
The second argument is optional, so, solve(Equate(A,B)); is enough.

simplify ... assuming ... should be enough.
Unfortunately an extra step is needed, for example:

arctan(cos(beta)/sin(beta));
simplify(convert(%, tan)) assuming beta>0, beta<Pi/2;

   Pi/2 - beta

Yes, it's a bug. With the col parameter the result is correct.

restart;
with(GraphTheory):
G1:=AddEdge(CycleGraph(4),{{1,3}}):
ChromaticNumber(G1);       # 2   wrong
ChromaticNumber(G1,'col'); # 3   correct
col;                       # [[2, 4], [3], [1]]

Edit. Strangely, all the methods give a wrong result:

for met in [hybrid, optimal, brelaz, dsatur, greedy, welshpowell, sat] do
  ChromaticNumber(G1, 'method'=met);
od;

Actually it seems that the bug is in GraphTheory:-GraphInfo:-GetAttrib
called by GraphTheory:-GetGraphAttribute.

For details and examples see  ?parallelism

P.S. It has happened to me many times that a Google search including a keyword "maple" is more successful than a search in Maple, even if the information is present in the help system.

I wrote a couple of years ago a procedure for this (multiple integrals over domains defined by polynomial inequalities).
https://www.maplesoft.com/applications/view.aspx?SID=154350

Actually, in that code the order of integration was chosen automatically.
Here is a version where the order is set by the user.

A more general but less efficient solution is

xx := (x,y,z) -> PUV(x,y,z)[1]:
yy := (x,y,z) -> PUV(x,y,z)[2]:
zz := (x,y,z) -> PUV(x,y,z)[3]:
fieldplot3d([xx,yy,zz], -1 .. 1, -1 .. 1, -1 .. 1, grid = [4, 4, 4], orientation = [160, 110, -100], arrows = `3-D`)

Of course here PUV is called 3 times instead of once. Maple should accept list- or vector-valued procedures [BTW, actually in fsolve too].

The exception indexing function is the only one which is not builtin. This seems to be the reason why imbricated exception tables do not work as expected.
It is possible to change its code, see:
showstat(`index/exception`);

but it could be tricky to find a working version.

Note that the sparse indexing works but of course unassigned entries will be 0 instead of an error.
In Maple 2019 (I see that unfortunately you use  Maple 2018) one can use e.g. the indexing sparse=undefined.
It is even possible to get errors for unassigned entries like this:

restart;
err:=proc() error "table exception" end:
b:=table(sparse='err()');
                 b := table(sparse = err(), [])
b[1]:=table(sparse='err()');
               b[1] := table(sparse = err(), [])
eval(b);
         table(sparse = err(), [1 = table(sparse = err(), [])])
b[1][2]:=12;
                         b[1][2] := 12
eval(b);
      table(sparse = err(), [1 = table(sparse = err(), [2 = 12])])
b[1][2];
                               12
b[1][3];
Error, (in err) table exception

You use a hackish (unorthodox) way to save random variables.
Why don't you save a proc (or module) to generate them?
Note that e.g. RandomVariable(Normal(0, 1)) returns simply a name _Rnn and _Rnn has attributes.
But the attributes are actually names for modules which contain the infos about the RV.
So these modules should be saved too.
But even so, I am not sure whether this will work, maybe Statistics does some other things too ...

This is a standard numeric behavior, covered in any numerical analysis textbook.
When a := 1+exp(-64) is computed with Digits=32, the accuracy for a-1 will be only 4 digits!  So the final result result will have only 4 correct digits.
Normally, Digits must be increased. But in this case we will be able to use Digits=32 and obtain max accuracy, with some maths. E.g.

convert(series(ln(1+x),x,3),polynom);
evalf[32](eval(%, x=exp(-64)));

        x - 1/2*x^2

        1.6038108905486378529760870340137*10^(-28)

The overlapping order is a bit more complicated. See here.

restart;
plots:-display(
seq(plot((x^2+x*y+y^2)*sqrt(((1-y^2)/(x+y))^2-(x-y)^2), x=0..1, color=HUE(y*0.66)), y=0..1,0.1),
legend=[seq(0..10)], labels=["x","z"]);

So, you want to hack VolumeOfRevolution:-ModuleApply, substituting 2*Pi with A and then animate. But 2*Pi appears in several places.

Why don't you try a standard animation with plot3d and animate?
 

For n >= 2, a(n) satisfies a(n)=2*a(n-1)+1,  [ a(n) - a(n-1) = 1 + a(n-1)], so:

rsolve({a(n)=2*a(n-1)+1, a(2)=2}, a);
seq(%, n=2..15);

       (3*2^n)/4 - 1

       2, 5, 11, 23, 47, 95, 191, 383, 767, 1535, 3071, 6143, 12287, 24575