I am trying to write a simulation in MATLAB for a $1D$ heat flow problem in which there are two materials in contact with each other but with different thermal properties. The results are unstable so I am trying to see what I'm doing wrong. Here is my problem $$ \dfrac{\partial T}{\partial t} = k_s\dfrac{\partial^2 T}{\partial x^2}, -1<x<0 $$ $$ \dfrac{\partial T}{\partial t} = k_i\dfrac{\partial^2 T}{\partial x^2}, 0<x<1 $$ with boundary conditions $$ T(0^-, t) = T(0^+, t) $$ $$ \lambda_sT_x(0^-, t) = \lambda_iT_x(0^+, t) $$ $$ T(-1, t) = 1, T(1, t) =0 $$
The numerical schemes I used for the interior points were normal implicit finite difference methods with $d_s = k_s \Delta t/\Delta x^2$ and $d_i = k_i \Delta t/\Delta x^2$:
$$ -d_sT_{j+1} ^{n+1} + (2d_s+1)T^{n+1}_{j} - -d_sT_{j-1}^{n+1} = T^n_j $$ for $j < 0$ and $$ -d_iT_{j+1} ^{n+1} + (2d_i+1)T^{n+1}_{j} - -d_iT_{j-1}^{n+1} = T^n_j $$
for $j > 0$. Whenever we have $j+1$ or $j-1$ equal to $0$, we use the average of $ d_s, d_i$ to represent the transfer coefficient at the connection point. To satisfy the continuity equation $T(0^-, t) = T(0^+, t)$. I just used the same point for both finite difference stencils. However, I tried implementing the continuity of flux boundary condition with the following logic
$$ \lambda_sT_x(0^-, t) = \lambda_iT_x(0^+, t) \approx \lambda_s \dfrac{T^{n+1}_{0} - T^{n+1}_{-1}}{\Delta x} = \lambda_i \dfrac{T^{n+1}_{1} - T^{n+1}_0}{\Delta x} \implies $$ $$ -\lambda_i T^{n+1}_1 + (\lambda_s + \lambda_i) T^{n+1}_0 - \lambda_s T^{n+1}_{-1} = 0 $$
However, in the case I considered where $\lambda_i = \lambda_s$, this scheme became unstable and exploded. I think my math is wrong about the continuity of flux, or just using the average of the two coefficients as the transfer at $0$. Does anyone have any insight on how to implement this?