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

@geri23 What I wrote was:

#######################################
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;
#The specific example:
params:={C1=25,C2=25,U12=40,U21=40,U13=20,U23=20,Ta=20,Tb=10,H1=(t->600), H2=(t->300),T3=(t->-5)};
#You can make H1, H2, T3 more interesting by making them not constant, e.g. making them periodic functions with period 24 hours. You need to think about what your unit of time is (second, minute, hour or what).
#The concrete system including initial conditions:
sys:=eval({eq1,eq2, ics},params);
#######################################
To make H1, H2, and T3 nonconstant functions just change params to reflect that. In params they are already defined as functions. Supposing e.g. that you want H1(t) to be 500+100*sin(t/T*2*Pi) for some concrete period T. Then you just replace H1=(t->600) by H1=(t->500+100*sin(t/T*2*Pi) ).

A numerical solution of an ODE is basically a solution solved in discrete time. So if your question is how to solve an ode numerically roughly speaking by hand, then look up the simplest possible method: Euler's method. You can even force Maple to use Euler's method when solving numerically by adding the optional argument method = classical[foreuler]. See ?dsolve,numeric,classical

A strongly related question was asked in

http://www.mapleprimes.com/questions/144357-Can-All-4x4-Matrices-Be-Written-As-A#comment144416