The computation of the analytical solution of the steady temperature distribution in multilayered media can become numerically unstable if there are different longitudinal (i.e., the directions parallel to the layers) boundary conditions for each layer. In this study, we develop a method to resolve these computational difficulties by approximating the temperatures at the junctions step-by-step and solving for the thermal field separately in only the single layers. First, we solve a two-layer medium problem and then show that multilayered media can be represented as a hierarchy of two-layered media; thus, the developed method is generalized to an arbitrary number of layers. To improve the computational efficiency and speed, we use varying weighting coefficients during the iterations, and we present a method to decompose the multilayered media into two-layered media. The developed method involves the steady-state solution of the diffusion equation, which is illustrated for 2D slabs using separation of variables (SOV). A numerical example of four layers is also included, and the results are compared to a numerical solution.