Use the Crank-Nicholson method to set up tridiagonal system of equations to
solve for ts at time n+1, where the temperature equation for layer i is
r_i = a_i [ts_i-1] n+1 + b_i [ts_i] n+1 + c_i [ts_i+1] n+1