solve({coeffs(expand(u-f), x)});

However, you have five variables and only four equations, so you won't get a complete solution.

In the following three function definitions

Ty:= y0-> T*diff((r-> rhs(r[-1]))@pds:-value(T(x, y), x = 0), [1], [y0]); 
Ty2:= y0-> T*diff((r-> rhs(r[-1]))@pds2:-value(T(x, y), x = 0), [1], [y0]); 
Ty3:= y0-> T*diff((r-> rhs(r[-1]))@pds3:-value(T(x, y), x = 0), [1], [y0]);

you need to use numeric differentiation, fdiff, rather than diff (which is only intended for symbolic differentiation). I don't have time to go into the details right now, but it looks like you were using fdiff syntax anyway.

You may also need to reduce the precision of the numeric differentiation for it to return results. Four digits of precision should be enough for any plot. I think that numpoints= 1000 is excessive given the number of mesh points used by pdsolve. You did say spacestep = .25. And I'm not saying that your mesh should be ~1000 points in the space dimension! That would probably take too long to solve.

Also, to ensure that you get multiplication, make sure that you enter "T space bar diff" or T*diff. In the last two of the above three cases, you actually had the single word Tdiff, which is meaningless, and I corrected it. But I'm fairly certain that your attempt to use T as a functional operator like that is incorrect. T and T(x,y) are just symbols; you need to use pds:-value again to access/create a numeric functional operator that corresponds to T.

You can either use

if Equal(Row(A,2), ZeroVector[row](6)) then ....

or (my preference)

if andmap(`=`, A[2, ..], 0) then ....

You can localize the value Digits to a particular call to evalf like this:

evalf[1500](Pi);

This avoids the need to explicitly set the value of Digits.

A lot of your computation can be avoided: If x is a float, then op(1, x) is its base-10 integer mantissa. For example, op(1, 3.14) is 314. Integers can be easily split into their digits with convert(..., base, 10). This works for any base if you just change the 10.

Putting it all together, your computation can be reduced to the single line:

Statistics:-Histogram(convert(op(1, evalf[1500](Pi)), base, 10), discrete, thickness= 9);

The appearance of the plot can be improved with some options:

Statistics:-Histogram(
   convert(op(1, evalf[1500](Pi)), base, 10),
   discrete, thickness= 9, view= [-1..10, 0..0.12], axis[1]= [tickmarks= [$0..9]]
);

The following array of three animations (one for each axis/angle) shows clearly a consistent and symmetric pattern, although it's perhaps idiosyncratic. The pattern is:

• Increasing values of the first angle specified in the orientation option causes clockwise rotation (when viewed from the positive end of the axis) about the third axis.
• Increasing values of the second angle causes clockwise rotation about the second axis.
• Increasing values of the third angle causes clockwise rotation about the first axis.

In each animation, the sole tick mark on each axis is at the extreme positive end and is the axis's standard name (x, y, or z). The positive end of the color bar along each axis is violet.

TorusAndAxes:= plots:-display(
   plots:-tubeplot(
      [cos(t),sin(t),0], t= -Pi..Pi, radius= 1/3, tubepoints= 8, color= t, numpoints= 16
   ),
   seq(
      plots:-tubeplot(
         axis, t= -1.7..1.7, radius= 1/8, color= HUE(.85*(t+1.7)/3.4), style= patchnogrid,
         tubepoints= 8, numpoints= 16
      ),
      axis= combinat:-permute([t,0,0])
   ),
   seq(axis[k]= [thickness= 2, tickmarks= [2=[x,y,z][k]]], k= 1..3)
):
   
plots:-display(
   Array(
      [seq(
         plots:-animate(
            plots:-display,
            [
               TorusAndAxes,
               orientation= ORI,   
               title= typeset(orientation= round~(ORI))
            ],
            d= 0..35, frames= 36, paraminfo= false
         ),
         ORI= combinat:-permute([10*d, 9, 9])
      )]
   ),
   axes= normal, thickness= 0,
   axesfont= [TIMES,BOLDITALIC,18],
   view= [(-2..2) $ 3],
   lightmodel= light1   
);

The animation is three copies of this thing, side by side:

The animation should be viewed in continuous mode at a frame rate of 1 per second.

Digits:= 50;
dsolve({sys2});
fsolve(rhs(%));

The value of Digits is just a guess. You just have to keep increasing it until the answer stabilizes.

It depends on what you want to see. If you want to see t on the horizontal axis vs. z on the vertical axis, then use the option scene= [t,z].

The command is 

A:= LinearAlgebra:-GenerateMatrix(convert(x, list), convert(u, list))[1];

The Equate command isn't necessary. Using convert(..., set) rather than convert(..., list) may cause the rows and/or columns of A to be returned in an unexpected order.

You can add the option orientation= [270, 0] to any 3-D plot command to get the view straight along the z-axis with the x- and y-axes in their standard 2-D orientation. If you're trying to understand the bounds of a double integral, this is the best view. 

What you want is not implicit differentiation in the ordinary sense of the term, there being no equation that implicitly defines alpha as a function of t. You just want alpha to be a function of t without the t appearing explicitly. The command for that is PDEtools:-declare.

Maple's equivalent of C++'s continue is next.

To do what you want requires a procedure (aka a function), and it can't be readily accessed through subscript notation:

a:= n-> n^2:
a(5), a(70);

If you really, really want to use a[5] instead of a(5), it is possible by creating an object for infinite sequences and overloading the `?[]` operator. This is probably not something for a new user to try.

Here's a very clean and robust way to do it:
 

Download UpdateStandings.mw

If x is between 0 and 1, then r(2x - 1) is between -r and r, and the transformation preserves uniformity. So the short Answer to your Question is that a correct procedure is

MyRand:= (M::posint, r::positive)-> ['r*(2*rand()/1e12 - 1)' $ M]:

But there's an easier way to use rand: When rand is passed a floating-point range, it returns a procedure that returns a random float in that range:

MyRandGenGen:= (r::positive)-> rand(-evalf(r)..evalf(r)):
MyRandGen:= MyRandGenGen(2):
#Test:
'MyRandGen()' $ 9;

1.48951387270106, -1.53779637227213, -.78650701082674, -1.93849247788427, -.50767257251262, 1.44464713288912, .80493082902846, 1.62849297809254, 1.18768745874873

If you want to generate long lists of floats, I think that the following is more efficient (but I haven't tested the efficiency):

MyRandGen:= (M::posint, r::positive)-> 
   RandomTools:-Generate(list(float(range= -r..r, method= uniform), M))
:
f:= MyRandGen(600, 2):

There are four dimensions to display: u, r, theta, and z. You've decided to handle u with a color range. A static (non-animated) plot can't show all the points in a three-dimensional space. We can get around this by showing an animation of a sequence of concentric cylindrical shells. So theta and z will be spatial dimensions, and r will be a time dimension.

restart:

#Example function
u:= (r, theta, z)-> r + theta + z:

#Example ranges:
Ranges:= Record(r= 0..2, theta= -Pi..Pi, z= -2..2, hue= 0..0.85):
#Use the hue range to tweak the color range. The colors are on the 
#standard visual spectrum: red = 0, violet= 1.  

#We need this module to collect the min and max of u so that the color
#scale is consistent through the whole animation.
U:= module()
export 
   umin:= infinity, umax:= -infinity,
   ModuleApply:= proc(r, theta, z)
   local ret:= u(r, theta, z);
      if not ret::realcons then return 'procname'(args) end if;
      if ret < umin then umin:= ret end if;
      if ret > umax then umax:= ret end if;
      ret
   end proc,
   Init:= proc() (umin,umax):= (1,-1)*~infinity end proc
;
end module:
 
#This first plot is simply to collect the min and max of u. 
#There's no need to see the plot.
U:-Init():
plots:-animate(plot3d, [U(r, theta, z), theta= Ranges:-theta, z= Ranges:-z], r= Ranges:-r):

plots:-animate(
   plot3d,
   [
      [r, theta, z], theta= Ranges:-theta, z= Ranges:-z, 
      color= op(2, Ranges:-hue) + `-`(op(Ranges:-hue))*(u(r,theta,z) - U:-umin)/(U:-umax - U:-umin),
      coords= cylindrical,
      #optional arguments:
      style= surface, axes= frame, labels= [``,``,'z'], grid= [36,40], lightmodel= light3,
      orientation= [45,60], thickness= 0, transparency= 0     
   ],
   r= Ranges:-r,
   #optional arguments:
   frames= 25
);

The images will be sharper and you'll be able to slow down the frame rate when you do this in a Maple worksheet. My color gradation is over the standard visual spectrum from red to blue. If you want to use strictly blends of pure red and pure blue, that can be done also. If your function u is strictly numeric (such as would be returned by pdsolve(..., numeric)), that will require a small adjustment to the above code.