**Background.** Finding analytical solutions to the problems of thermal conductivity with variable physical properties of the medium by classical analytical methods is very complicated mathematically. The known expressions representing complex infinite series including two types of Bessel functions and gamma-functions are, in fact, numerical as they require a numerical solution to complex transcendental equations with eigenvalues of the boundary problem. Such solutions can hardly be used in engineering applications, especially in cases when a solution to a certain problem is only an intermediate stage in other problems (such as thermoelasticity and control problems, inverse problems, etc.) which can be solved effectively only by finding analytical solutions to the initial problems. Therefore, an urgent problem now is to develop new methods of obtaining analytical solutions to the abovementioned problems, at least approximate ones.

**Materials and methods.** The study employed methods of additional boundary conditions and additional unknown functions in the integral method of heat balance.

**Results.** High-precision approximate analytical solutions to the transient heat conduction problem with nonhomogeneous physical properties of the medium for an infinite plate under symmetric boundary conditions of the first type have been obtained. The initial problem for partial differential equations is reduced to two problems in which ordinary differential equations are integrated. Additional boundary conditions are defined in such a way that their fulfillment in accordance with the new method is equivalent to the result of solving the initial partial differential equation at the boundary points and at the temperature perturbation front (for the first stage of the process).

**Conclusions.** By combining methods with finite and infinite heat propagation rate we have been able to obtain high-precision analytical solutions for the whole time range of the unsteady process including its small and ultra small values. The solutions look like simple algebraical polynomials not including special functions (Bessel, Legendre, gamma-functions and others). Since it is not necessary to directly integrate the initial equations by the space variable and to reduce them to ordinary differential equations with additional unknown functions, the considered method can be used for solving complex boundary problems in which differential equations do not allow distinguishing between the variables (into nonlinear, with linear boundary conditions and heat sources, etc.).