Now what went before causing you to write "so i got this"?

How do you make sense out of that question, or is it one? "So i got rhis ... "

@LEETZ You can use value:

uval:=subs(s:-value(output=listprocedure),u(r,t));
#The first command gives the same as s:-animate:
plots:-animate(plot,[uval(r,t),r=0..1],t=0..5);
#Without using cylindrical coordinates:
plots:-animate(plot3d,[uval(sqrt(x^2+y^2),t),x=-1..1,y=-sqrt(1-x^2)..sqrt(1-x^2)],t=0..5);
#Using cylindrical coordinates:
addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
plots:-animate(plot3d,[uval(r,t),r=0..1,theta=0..2*Pi,coords=z_cylindrical],t=0..5);

@LEETZ You can use value:

uval:=subs(s:-value(output=listprocedure),u(r,t));
#The first command gives the same as s:-animate:
plots:-animate(plot,[uval(r,t),r=0..1],t=0..5);
#Without using cylindrical coordinates:
plots:-animate(plot3d,[uval(sqrt(x^2+y^2),t),x=-1..1,y=-sqrt(1-x^2)..sqrt(1-x^2)],t=0..5);
#Using cylindrical coordinates:
addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
plots:-animate(plot3d,[uval(r,t),r=0..1,theta=0..2*Pi,coords=z_cylindrical],t=0..5);

@geri23 You may want to take a look at the help page:

?dsolve,events

Here is a simple version of a  thermostat model inspired by Joe Riel a few years ago:

dsys := {diff(T(t),t) = piecewise(b(t)=0, -k*(T(t)-Tu), -k*(T(t)-Tu)+F)
          , T(0) = T0
          , b(0) = 1
        };

L := dsolve(dsys, numeric
                , 'events' = [[T(t)=T1, b(t)=1],[T(t)=T2, b(t)=0]]
                , 'discrete_variables' = [ b(t) :: boolean ]
                , parameters=[Tu,T1,T2,T0,k,F]):
L(parameters=[10,18,22,7,1,13]);
plots:-odeplot(L, [[t,T(t)]], 0..10);

@geri23 If that happens then it likely illustrates that discrete methods may be unstable under certain conditions, in particular if the stepsize is too large.
Why do you insist on a fixed stepsize method. What is the point? Is it perhaps to exhibit the instability?

@geri23 Apparently your time unit is second (T is 24 hours). So you may want to solve the system using Euler's method with stepsize = 5*60 (5 minutes).

So you don't insist that the matrix should be a rotation matrix, but instead that the sum of the rows and the first two columns all be 1 and the determinant be 1?

with(LinearAlgebra):
R:=Matrix([[R11,R12,fx],[R21,R22,fy],[R31,R32,fz]]);
Determinant(R);
Transpose(R).R=IdentityMatrix(3); #You don't require this, or do you?

Why the second root? The ordering of the output may not be what you want:

for alpha from .4 by .1 to 1 do
   argument~([solve(x^4-alpha, x)]);
end do;

Number 2 root has argument Pi/2 in all but the last case, where it is Pi.

You could sort the output e.g. in one of these ways:
for alpha from .4 by .1 to 1 do
   sort([solve(x^4-alpha, x)],(x,y)->evalb(Re(x)<Re(y) or Re(x)=Re(y) and Im(x)<Im(y)));
end do;
#or
for alpha from .4 by .1 to 1 do
   sort([solve(x^4-alpha, x)],(x,y)->evalb(argument(x)<argument(y)));
end do;

The advantage of Markiyan's version is that the output is the sequence sought. A print statement just prints to the screen. You have to copy the result from the screen afterwards.
Using a loop the same could be done with a small change:

S:=NULL:
for alpha from .4 by .1 to 5 do
   S := S,solve(x^4-alpha, x)[2];
end do:

S;

The advantage of Markiyan's version is that the output is the sequence sought. A print statement just prints to the screen. You have to copy the result from the screen afterwards.
Using a loop the same could be done with a small change:

S:=NULL:
for alpha from .4 by .1 to 5 do
   S := S,solve(x^4-alpha, x)[2];
end do:

S;

@LEETZ 
You could do
s:-plot3d(t=0..10, r=0..1);

@LEETZ 
You could do
s:-plot3d(t=0..10, r=0..1);

@Xucu Maybe something like this:
restart;
eq:=diff(y(t),t,t)+u(t)-y(t)=0;
res:=dsolve({eq,y(0)=0,D(y)(0)=1,u(0)=0,n(0)=1},numeric,output=listprocedure,
  discrete_variables=[u(t)::float,n(t)::integer],
  events=[[t=[0,5],[u(t)=y(t)+diff(y(t),t)+(-1)^n(t)*0.02,n(t)=n(t)+1]]]);

plots:-odeplot(res,[[t,y(t)],[t,u(t)],[t,n(t)]],0..50);

@Xucu Maybe something like this:
restart;
eq:=diff(y(t),t,t)+u(t)-y(t)=0;
res:=dsolve({eq,y(0)=0,D(y)(0)=1,u(0)=0,n(0)=1},numeric,output=listprocedure,
  discrete_variables=[u(t)::float,n(t)::integer],
  events=[[t=[0,5],[u(t)=y(t)+diff(y(t),t)+(-1)^n(t)*0.02,n(t)=n(t)+1]]]);

plots:-odeplot(res,[[t,y(t)],[t,u(t)],[t,n(t)]],0..50);