The objective of the proposed research is to address the multi-physics design challenges associated with the advanced packaging and heterogeneous integration of integrated voltage regulators (IVRs). IVRs, which provide high power delivery performance, are becoming increasingly common in modern computing systems. However, their advantages come with significant multi-physics design complexities. In particular, high-voltage conversion IVRs typically require the integration of passive components, such as inductors, into the package. For the IVR to function effectively, these package-embedded inductors must meet several design criteria, including appropriate inductance, thermal resistance to ambient temperature, saturation current, and compact size. However, there is often a trade-off between inductance and saturation current—package-embedded inductors tend to either have high inductance with low saturation current or vice versa. Furthermore, these inductors are subject to Joule heating and magnetic core losses, resulting in high thermal resistance and elevated temperatures, which can compromise the efficiency of the IVR. To tackle these challenges, computational modeling and simulation are essential. However, traditional time-domain modeling methods face significant computational complexity due to their repetitive time-stepping nature and the involvement of multiple physical domains. To overcome both the issue of low saturation current in embedded inductors and the complexities of multi-physics modeling, we propose two innovative solutions: dual-core coupled spirals and an augmented computational method. The dual-core coupled spiral topology features two spiral coils sandwiched between magnetic cores, physically separated by an air gap. This gap increases the saturation current capacity of the inductor, improving its performance. The augmented computational method, on the other hand, employs a frequency-domain approach to multi-physics modeling, enabling a unified tool for both transient time-domain electrical simulations and seamless steady-state thermo-mechanical analysis. This method uses Fourier series[1]based Toeplitz matrices and vectors to represent the time-varying nature of elements and signals within the system. Additionally, the complexity of the augmented computational method can be further reduced through parallelization on CPU/GPU architectures, significantly improving efficiency and enabling faster simulations.