@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

@trace I would begin by making m concrete. As many examples as possible.

@trace I would begin by making m concrete. As many examples as possible.

@trace Near the bottom of p.18 the authors make a change of variable. So that has to be done.
Here are the beginnings for the case considered m(r) = n-1.
I start as before.
restart;
f := (x, y, z)-> -1-z-3*y+x^2-x*z;
g := (x, y, z)-> x*z/m(z/y)-2*z*y+4*y+x*y;
h := (x, y, z)-> -x*z/m(z/y)-2*z^2+4*z;
m:=s->n-1;
res := solve({f(x, y, z) = 0, g(x, y, z) = 0, h(x, y, z) = 0}, {x, y, z});
pts:=map(subs,[res],[x,y,z]);
J := unapply(VectorCalculus:-Jacobian([f(x, y, z), g(x, y, z), h(x, y, z)], [x, y, z]), x, y, z):
use ev=LinearAlgebra:-Eigenvalues in
  for p in pts do
  ev(J(op(p)))
  end do
end use;
#For convenience I'm using x1,x2,x3 instead of x,y,z:
sys := diff(x[1](t), t) = f(x[1](t),x[2](t), x[3](t)), diff(x[2](t), t) = g(x[1](t),x[2](t), x[3](t)), diff(x[3](t), t) = h(x[1](t),x[2](t), x[3](t));
p:='p':
alias(P=add(p[i](t)^2,i=1..3)); #To make output shorter
varchange:={seq(x[j](t)=p[j](t)/(1-sqrt(add(p[i](t)^2,i=1..3))),j=1..3)};
PDEtools:-dchange(varchange,{sys},[seq(p[i](t),i=1..3)]);
solve(%,{seq(diff(p[i](t),t),i=1..3)});
psys:=collect(%,[p[1](t),p[2](t),p[3](t)],distributed,factor);
#This is the system in Poincare coordinates. You need to specify n, solve numerically and show plots in the (p1,p2)-plane.

@Markiyan Hirnyk The answer returned by fsolve doesn't mean anything here. Try

evalf[40](eval(J,{z=0.1e-2,x=-12.30640295}));
evalf[40](eval(J,{z=0.1e-2,x=-12.4}));

The problem surely is that J itself is extremely small at those x-values. The "roots" are not roots at all.

My argument was and remains: Since diff(J,x) < 0 for x<0 and limit(J,x=-infinity) = 0 the equation J = 0 has no negative root for x for any given value of z. 

@Markiyan Hirnyk The answer returned by fsolve doesn't mean anything here. Try

evalf[40](eval(J,{z=0.1e-2,x=-12.30640295}));
evalf[40](eval(J,{z=0.1e-2,x=-12.4}));

The problem surely is that J itself is extremely small at those x-values. The "roots" are not roots at all.

My argument was and remains: Since diff(J,x) < 0 for x<0 and limit(J,x=-infinity) = 0 the equation J = 0 has no negative root for x for any given value of z. 

@orawajar You can use the Routh-Hurwitz criterion. For a polynomial with real coefficients this gives necessary and sufficient conditions for all roots to have negative real parts.
In case of a polynomial of 4th degree the conditions are:

p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
for i from 0 to 4 do a[i]:=coeff(p,lambda,i) end do;
All the coefficients have to be positive and furthermore
a[3]*a[2]>a[4]*a[1];
a[3]*a[2]*a[1]>a[4]*a[1]^2+a[3]^2*a[0];
#Not surprising then that the following returns FAIL:
PolynomialTools:-Hurwitz(p,lambda);

@orawajar You can use the Routh-Hurwitz criterion. For a polynomial with real coefficients this gives necessary and sufficient conditions for all roots to have negative real parts.
In case of a polynomial of 4th degree the conditions are:

p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
for i from 0 to 4 do a[i]:=coeff(p,lambda,i) end do;
All the coefficients have to be positive and furthermore
a[3]*a[2]>a[4]*a[1];
a[3]*a[2]*a[1]>a[4]*a[1]^2+a[3]^2*a[0];
#Not surprising then that the following returns FAIL:
PolynomialTools:-Hurwitz(p,lambda);

@Axel Vogt Right, if the normal/expanded line wasn't there things work out fine in this case.
However, with that line:

eval(1/(2^s-1), s = x - 2*Pi*I*3/log(2));
normal(%,expanded);
series(%, x=0);
coeff(%, x, -1);
Result 0.

@Alejandro Jakubi I forgot to take care of the case s = infinity. Translating to zero makes the code somewhat simpler, but not much. Set x,x0:=op(a) and after assigning to g in the infinity case set x0:=0, then make the trivial adjustments.