We have two heat reservoirs 0 and 1, a channel between them of conductance P, and another channel of conductance Q between reservoir 1 and constant temperature 0.
(Conductance is the reciprocal of resistance.)
The rate of heat flow for the first channel is P times the temperature difference T_{2} − T_{1}.
The heat capacity of the reservoir j is C_{j}.
The total heat in reservoir j is T_{j}C_{j}.

CWhen we get to actually computing we will assume that P = 2, Q = 1, C_{0}dT_{0}/dt = P(T_{1}− T_{0}) C_{1}dT_{1}/dt = P(T_{0}− T_{1}) − QT_{1}

This falls into the more general form

C_{i}dT_{i}/dt = Σ_{j}c_{ij}T_{j},
where c_{00} = −P; c_{01} = P;
c_{10} = P; c_{11} = −P−Q.

and the naïve code to compute the answer.

Letting M_{ij} = c_{ij}/C_{i}
we can rewrite:

dT_{i}/dt = Σ_{j}M_{ij}T_{j}

Rewriting in Matrix form we get

dT/dt = MTHere is a curious scheme that serves for some purposes but is not what we want to study here.

Having invoked matrix notions, we go a step further and invoke vector spaces.
We view our set of temperatures T as a vector whose components are the temperatures of the individual components,
Our equation tells us how T changes with time, depending on T itself.
If we found another coördinate system in this space in which the matrix M' for our transformation was diagonal then in those coördinates we would have n independent equations each of the form
dT'_{i}/dt = D_{i}T'_{i} where T'_{i} and D_{i} are real numbers and D_{i} is the ith element on the diagonal of M'.
Finding such a transformation is the classic eigenvector problem and the D_{i}’s are the eigenvalues.
The matrix for the new coördinate system is M' = G^{−1}MG where G is a matrix whose columns are the eigenvectors.

Our new coördinate system uses the eigenvectors as a basis; F_{j} is the jth basis vector of the new coördinate system and we can write any solution to these new differential equations as Σ_{j}a_{j}F_{j}e^{Djt} for some appropriate set of a_{j}’s.
This Algol 68 program finds G and prints G^{−1}MG: a diagonal matrix.
The program id “g” denotes G which names a row of column vectors each of which is an eigenvector.
We speak of a *new coördinate system* whose basis vectors are these eigenvectors.
Choose an eigenvector x; expressed in the new coördinates it is <0, … 0, 1, 0, … 0>.
Gx names the same vector in the old coördinates.
MGx names the vector (λx) stretched in the same direction, expressed in the old coördinates.
G^{−1}MGx names the stretched vector in the new coördinates.
G^{−1}MG should thus be diagonal!
That is indeed the last matrix printed and it is diagonal with the eigenvalues on the diagonal!

After printing G^{−1}MG the values Σ_{j}a_{j}F_{j}e^{Djt} are printed in the old coördinates, for successive times.
The a_{j} are chosen to fit the initial conditions.
The answers agree with the simpler calculations.

This C code captures the business end of the above. It omits values calculated only to verify the math. It relies on the routine slv here and routine eigen here.

The line below `main()` has all the inputs.

In the equation C_{i}dT_{i}/dt = Σ_{j}c_{ij}T_{j},
“c” is `c` of the program.
The `d`’s of the program are the reciprocals of the “C”’s in the equation.
The definition of n must also be adjusted.
`sl[j][n]` is the a_{j} in the description of the logic above.

Here is an alternate derivation of some of the logic.