Are you sure that the original had a z-axis label? A basic plot3d command such as

plot3d(x^2+y^2, x= -2..2, y= -2..2);

doesn't produce a z-axis label unless you explicitly include one via, for example,

plot3d(x^2+y^2, x= -2..2, y= -2..2, labels= ['x', 'y', 'z']);

If that doesn't cure your problem, please upload the code that generates the plot.

fsolve returns a sequence of roots (when solving a polynomial) rather than a set of roots, but select expects a set. So, replace fsolve(r, x[2]) with {fsolve(r, x[2])}. When you do that, you're just going to get the empty set on each loop iteration because r has no roots between 0 and 1; indeed, it has no real roots at all. To see the complex roots, use fsolve(r, x[2], complex) without the select.

The code in the attached file runs without error, so I don't know what to correct. Perhaps you simply mean that you want to see the results? For that use

eval(Result_Eva_GB_X3_X2);

after the loop.

The results of code inside loops is ordinarily not shown. That can be adjusted by setting the value of printlevel (see ?printlevel). One needs to be careful with this because it's very easy for the amount of information displayed to become overwhelming, and once it has started displaying there's usually no way to stop it other than killing the whole kernel.

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.