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

@geri23
1. Yes to both questions.
2. Comparison is easiest if you choose output=listprocedure in dsolve/numeric rather than the default (procedurelist).
Here is the previous code with some additions. (Edited May 31)

restart;
eq1:=C1*diff(T1(t),t)=U12*(T2(t)-T1(t))+U13*(T3(t)-T1(t))+H1(t);
eq2:=C2*diff(T2(t),t)=U21*(T1(t)-T2(t))+U23*(T3(t)-T2(t))+H2(t);
ics:=T1(0)=Ta, T2(0)=Tb;
T:=15: #Choosing period T = 15 for H1:
params:={C1=25,C2=25,U12=40,U21=40,U13=20,U23=20,Ta=20,Tb=10,H1=(t->600+100*sin(2*Pi/T*t)), H2=(t->300),T3=(t->-5)};
sys:=eval({eq1,eq2, ics},params);
res:=dsolve(sys);
sol:=subs(res,[T1(t),T2(t)]);
plot(sol,t=0..100,caption="Analytical solutions",legend=["T1","T2"],legendstyle=[location=right]);
resnum:=dsolve(sys,numeric); #Uses the default numerical method rkf45
plots:-odeplot(resnum,[[t,T1(t)],[t,T2(t)]],0..100,caption="Numerical solution by the RKF45-method",legend=["T1","T2"],legendstyle=[location=right]);
plots:-odeplot(resnum,[[t,T1(t)-sol[1]],[t,T2(t)-sol[2]]],0..10,legend=["RKF45 - analytic for T1","RKF45 - analytic for T2"],legendstyle=[location=bottom]);
resEuler:=dsolve(sys,numeric,method=classical[foreuler],stepsize=.1); #Euler's method
plots:-odeplot(resEuler,[[t,T1(t)],[t,T2(t)]],0..100,caption="Numerical solution by Euler's method.\nStepsize=0.1",legend=["T1","T2"],legendstyle=[location=right]);
plots:-odeplot(resEuler,[[t,T1(t)-sol[1]],[t,T2(t)-sol[2]]],0..10,legend=["Euler - analytic for T1","Euler - analytic for T2"]);
#Now using output=listprocedure:
resEulerLP:=dsolve(sys,numeric,method=classical[foreuler],stepsize=.1,output=listprocedure);
#Extracting the numerical procedures for T1 and T2 obtained by Euler's method with stepsize=0.1:
TE1,TE2:=op(subs(resEulerLP,[T1(t),T2(t)]));
plot([TE1(t),sol[1]],t=0..10,caption="Analytical solution and Euler solution for T1",legend=["Euler","Analytical"],legendstyle=[location=right]);