Phase-change materials (PCMs) excel in storing significant thermal energy through the latent heat of fusion during phase changes. However, they often suffer from low thermal conductivity, requiring the incorporation of materials with higher thermal conductivity to overcome this limitation. In this study, we employ the level-set method to topologically optimize the distribution of both PCM and high thermal conductivity materials within heat sinks.The heat transfer process is modeled as an unsteady diffusion problem, with PCM integrated using the apparent heat capacity method. The finite element method, implemented using FEniCS, is utilized to compute the temperature distribution at each time step. The optimization problem is solved using ParaLeSTO, a modular topology optimization software written in C++, seamlessly interfacing with FEniCS through its Python interface.To compute boundary point sensitivities for level-set topology optimization, the discrete adjoint method is employed, combining a perturbation scheme with automatic differentiation. The proposed methodology is first applied to a two-dimensional example of temperature oscillation minimization with a heat flux boundary condition, drawn from existing literature. Subsequently, a two-dimensional example with multiple volumetric heat sources is studied, followed by an extension to a three-dimensional design domain.Comparisons with optimized designs without PCM, serving as a reference, highlight the superior thermal performance of heat sink designs with PCM. The results showcase a maximum temperature reduction of up to 42%, a mean temperature reduction of up to 42%, and a temperature variance reduction of up to 94%. This research contributes to heat sink design by offering discrete solutions that significantly improve thermal performance as compared to designs without PCM. Its impact extends to addressing crucial challenges in thermal management across various applications, including automotive cooling systems, aerospace components, and consumer electronics.