@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);

You made a silly mistake (nice to see that I'm not the only one!).

You made a silly mistake (nice to see that I'm not the only one!).

@mehdi_mech Use output=listprocedure:
Here illustrated without scaling (per our former discussion) and with your original equations, which only included a[i,j] for i,j = 1..3:
I assign to A[i,j] (i,j=1..3) the numerical procedures for computing a[i,j](t).

MMM := dsolve(AA,numeric,maxfun=10^7,output=listprocedure);
Alist:=[seq(seq(A[i,j],i=1..3),j=1..3)];
alist:=[seq(seq(a[i,j](t),i=1..3),j=1..3)];
assign(Alist=~subs(MMM,alist));
Alist;
eval(A[1,1]);
plots:-odeplot(MMM,[t,a[1,1](t)],0..1e-4);
u:=A[1,1](s)*x+A[2,1](s)*x^2+A[3,1](s)*x^3;
plot ( eval(u,s=0.0004),x=0..1);
Comment:
Since there are so many procedures I'm doing it in a somewhat complicated way. Had there only been 2 unknowns a1(t) and a2(t) I would have done like this:
A1,A2:=op(subs(MMM,[a1(t),a2(t)]));

@mehdi_mech Use output=listprocedure:
Here illustrated without scaling (per our former discussion) and with your original equations, which only included a[i,j] for i,j = 1..3:
I assign to A[i,j] (i,j=1..3) the numerical procedures for computing a[i,j](t).

MMM := dsolve(AA,numeric,maxfun=10^7,output=listprocedure);
Alist:=[seq(seq(A[i,j],i=1..3),j=1..3)];
alist:=[seq(seq(a[i,j](t),i=1..3),j=1..3)];
assign(Alist=~subs(MMM,alist));
Alist;
eval(A[1,1]);
plots:-odeplot(MMM,[t,a[1,1](t)],0..1e-4);
u:=A[1,1](s)*x+A[2,1](s)*x^2+A[3,1](s)*x^3;
plot ( eval(u,s=0.0004),x=0..1);
Comment:
Since there are so many procedures I'm doing it in a somewhat complicated way. Had there only been 2 unknowns a1(t) and a2(t) I would have done like this:
A1,A2:=op(subs(MMM,[a1(t),a2(t)]));