I have a perfectly working 2-D finite-difference thermal solver, using an alternating-direction implicit scheme (painfully written in Fortran 90 at my adviser’s insistence), that I have recently extended to 3-D. It is essentially a three step approximation to the Crank-Nicolson equation, which is shown to be unconditionally stable in homogeneous media with constant boundary conditions in the old text book in which I found it. Who the hell writes a thermal solver to simulate a single material with constant temperatures at the boundaries? Anyway, I have adapted the method to an inhomogeneous material structure, and, regardless of boundary conditions, I have exponential error growth (in magnitude...the actual error oscillates), which starts at the material interfaces. For almost prohibitively low time-steps, though still much larger than what would be required for an explicit method, the thermal solver actually appears stable and produces the correct results as far as I can tell, so I don't think there is any simple error in my implementation. It appears to me that the material inhomogeneities have broken the unconditional stability of the 3-D alternating direction implicit scheme, and I am having a very hard time finding relevant resources. Posting this here is probably a bit of a fishing expedition, but if this sounds familiar to anyone else, please elaborate on your experiences. Thanks in advance.