thermodynamics - How do I set up the tridiagonal matrix for a heat diffusion with layers of different thermal diffusivity?

I have Scala code that recreates the Crank-Nicolson solutions for the diffusion equations, and matches 'Excel for Scientists and Engineers' (Joe Billo, Wiley).

However, I would like to be able to replicate the results from H. Asan's 'Numerical computation of time lags and decrement factors' (including the decrement factor with very thin materials) and also the multi-layer positioning results from the BS2013 paper 'Optimizing insulation-thermal mass layer distribution from maximum time lag and minimum decrement factor point of view'.

Both papers suggest that Crank-Nicolson is at the heart of the approach, however there is a complexity with the treatment of energy flow at the outer and inner surface, and with different material layers, that are not covered in the simple diffusion treatments that I found.

How do I set up the matrix - or how do I connect a series of such matrices together?

I had thought that perhaps a simple treatment would give a given wall thickness a behavior like a delay factor and an amplitude scaling - but it does not seem to be possible to do this and sum the phase delays and multiply the decrements, because the results above suggest that the ordering is not transitive.

How can I proceed - given that I'm a lot better at coding than solving PDEs.