Thermal simulation of a Printed Circuit Board (PCB) can help identify potential overheating risks in the circuit. The proposed modeling method combines analytical temperature solutions and numerical approximations. Only Fourier-series analytical solutions related to the prepreg-layer surfaces need to be calculated, rather than the entire structure. Heat transfer through the lateral sides of a PCB is approximately considered as part of the compensated heat flux of the insulating-layer surface boundaries. Heat diffusion within or between metal layers is numerically approximated using the finite volume method. The core layer is treated as "thermally-thick". Temperature-dependent boundary conditions are considered through iterations. A test solver was developed based on the method. The modeling accuracy was validated by comparison with COMSOL Multiphysics for a four-layer structure with a moderate degree of discretization. Additionally, a PCB for generating DC 3.3V was designed, tested, and modeled, with the modeling results confirmed by the thermal images. The electro-thermal analysis of the distribution of electric potential and Joule heating in traces and vias was integrated into the PCB model. The layout maps of the PCB were further adjusted to reduce Joule heating in the output circuit, and the improvement on reducing the IR drop and hotspot temperature was examined.