Add eigenvalue decomposition to thermal modelling

[bcych]

• Get eigenvalue decomposition for Q_matrix * dt working.
• Divide through by fac = mean(abs(Q_matrix))
• Multiply by exp(Q_matrix * dt) * exp(fac)
• Check matrix is possible to factorize!

This should greatly improve numerical stability at large dt.

[bcych]

Set to low priority as multiprecision is now working properly as of commit 610688b