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

A positive hardware float (hfloat)  is "normally" in the range  10^(-308) .. 10^(308) [approx].
But  compromising precision, a subnormal representation (having the mantissa starting with zeros) allows smaller values up to about 5*10^(-324).
Your number 9.395*10^(-315) requires such a subnormal hfloat.
Note that evalhf  does conversions float -> hfloat -> float.
It seems that there is a bug in these conversions for subnormal hfloats. The bug seems to be very old.

restart;
interface(prettyprint=0):
Digits:=20:
evalhf(9.87654321*10^(-315));
-0.85702008631607784918e-314
evalhf( proc()9.87654321*10^(-315)end ());
0.
9.87654321*10^(-315); evalhf(%);
0.98765432100000000000e-314
-0.85702008631607784918e-314
k:=315: evalhf(9.87654321*10^(-k));
0.98765431971986665421e-314
for n from 305 to 330 do
  evalhf(9.87654321*10^(-n));
od;
0.987654321000016695e-304
0.987654321000016829e-305
0.987654321000016871e-306
0.987654321000017055e-307
0.98765432100001643399e-308
0.98765432099999568324e-309
0.98765432099994133602e-310
0.98765432099850854564e-311
0.98765432101431864631e-312
0.98765432095503076881e-313
0.98765431971986665421e-314
0.98765430588602857065e-315
0.98765456774082086651e-316
0.98765303613731875865e-317
0.98764216669311025122e-318
0.98763722603665183876e-319
0.98566096345328685259e-320
0.97824997876566815445e-321
0.98813129168249308530e-322
0.
0.
0.
0.
0.
0.
0.

At each step except the first, you take tau12, tau21 from the previous step.

Numeric only (of course).

ode:=((D@@2)(y))(r)+(D(y))(r)/r-sin(y(r))*cos(y(r))/r^2+sin(y(r))^2/r-sin(y(r))-sin(2*y(r)) = 0:
ds:=dsolve([ode, y(2)=1, D(y)(2)=1/2], numeric):
plots:-odeplot(ds, r=1e-6..100, color=red);

Nothing unusual. The integrator is very complex, so, equivalent functions could be integrated using different methods.

But it's easy to have the same unevaluated result using:
int(simplify(integrand2), x);  ##  :-)

In the absence of a better solution (the generated LaTeX code is not easy to manipulate), I'd replace in the worksheet e.g.

sys_set := convert(sys_vec, set);

with

sys_set := convert(sys_vec, set):
for eq in sys_set do print(eq) od; # display elements

You don't need Maple for this: z / |z|^2 = z is holomorphic, so the answer is 0.

The limit of a function F : A --> B, at a point u is defined when u is an accumulation point (= limit point = cluster point) of A.
In our case the standard choises for A, B are B = R and A = R \ {0} \ {1/(n*Pi) : n in Z\{0}}.
0 is an accumulation point for A, so the limit problem makes sense. It exists and equals 1 (using the standard limit : sin(t)/t, for t --> 0).

[A more techical note: the topology in A is by default the subspace topology induced by the natural topology in R; your neighbourhood refers to this topology].