Finding eigenvectors is fast, easy and efficient when the matrix is symmetric, but alas M is not symmetric unless the C_i are all the same. To get around this we note that c_ij is symmetric and adopt t_i = T_i√C_i as a temporary surrogates for the temperatures.
dT_i/dt = Σ_jM_ijT_j
(dt_i/dt)/√C_i = Σ_j(M_ijt_j/√C_j)
dt_i/dt = Σ_j((M_ij√C_i)t_j/√C_j)
dt_i/dt = Σ_j((M_ij√C_i)/√C_j)t_j
dt_i/dt = Σ_j(((c_ij/C_i)√C_i)/√C_j)t_j
dt_i/dt = Σ_j(c_ij/√(C_iC_j))t_j
dt_i/dt = Σ_jm_ijt_j where m_ij = c_ij/√(C_iC_j)

Note that the new matrix m is symmetric! c_ij is symmetric indeed for heat flow problems and most other physics problems where reciprocity reigns.

Let D be a diagonal matrix where D_ii = 1/C_i. We may now write the matrix equation m = DC (this is a different m!! Fix! ....) and use this code to compute eigenvectors and values for m.

We find g so that gmg⁻¹ is diagonal with diagonal elements d_j. If E_j is the jth basis vector in the new coordinates then our general solution to the equations is:
Σ_ja_jE_je^d_jt. The n reals a_j are the parameters to the solution.