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

@geri23 I have edited my code above to include captions and legends in order to make it somewhat more understandable. I hope that that will answer your questions.

"...because Maple won't accept diff(..., x$0)."
Carl, you can use a list as a second argument to diff:
diff(f(x),[]);
diff(f(x),[x$0]);
diff(f(x),[x$2]);

"...because Maple won't accept diff(..., x$0)."
Carl, you can use a list as a second argument to diff:
diff(f(x),[]);
diff(f(x),[x$0]);
diff(f(x),[x$2]);