  A class of models of heat transfer processes in a multilayer domain is
considered. The governing equation is a nonlinear heat-transfer equation with
different temperature-dependent densities and thermal coefficients in each
layer. Homogeneous Neumann boundary conditions and ideal contact ones are
applied. A finite difference scheme on a special uneven mesh with a
second-order approximation in the case of a piecewise constant spatial step is
built. This discretization leads to a pentadiagonal system of linear equations
(SLEs) with a matrix which is neither diagonally dominant, nor positive
definite. Two different methods for solving such a SLE are developed --
diagonal dominantization and symbolic algorithms.
