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.

If you use exact arithmetic (that means no decimal points), you'll see that the function ratio doesn't depend on its third parameter, l (ell). For example,

ratio(1, 2, L, 3, 4);

So the problem is not with the plot command. There's some mathematical reason why l cancels in your expression.

(For the sake of readability, it's probably not good to use lowercase l as a variable.)

Let's consider this block of your code:

if j-1>1 then
   if N <> Configz1[w](o, j+1) then DeltaBond:= DeltaBond+1
      if N = Configz1[w](o, j+1) then DeltaBond:= DeltaBond-1
else DeltaBond:= DeltaBond
end if:
end if:
end if:

Now let's look at it again, execpt that I'll use a logical indentation structure:

if j-1>1 then
    if N <> Configz1[w](o, j+1) then 
       DeltaBond:= DeltaBond+1
       if N = Configz1[w](o, j+1) then 
         DeltaBond:= DeltaBond-1
      else
         DeltaBond:= DeltaBond
      end if:
   end if:
end if:

Now do you see the problem? The third line needs a semicolon or colon at the end.

For BVPs, the equivalent of maxfun is maxmesh.

When or and seq or map are used togther, they can be replaced by a single command. In this case that's

ormap(`=`, R[..-2], R[i]);

This is more efficient than Kitonum's command because it aborts the processing as soon as the first true result is found.

There's also an andmap.