Objectives: Differentiated thyroid cancer (DTC) is usually treated with surgery and radioactive iodine therapy. However, in advanced stages, molecular changes can lead to a progressive loss of sensitivity to iodine, making the cancer resistant to radioactive iodine treatment (known as radioactive iodine-refractory, or RAIR). For cases of RAIR-DTC, alternative treatments such as tyrosine kinase inhibitors and immunotherapy are explored. To better understand the dynamics of these alternative treatments, a new mathematical model has been developed. In this model, Sorafenib represents the tyrosine kinase inhibitor, while Pembrolizumab is the immunotherapy agent under investigation. Methods: Two mathematical models were developed based on systems of ordinary differential equations to study differentiated thyroid cancer treated with (1) Sorafenib alone and (2) a combination of a Sorafenib and Pembrolizumab. The first model included three variables: drug concentration of Sorafenib, number of cancer cells, and number of immune cells with fourteen parameters. The second model additionally incorporated the concentration of the Pembrolizumab. Parameter values were estimated through simulations utilizing the fmincon optimization technique in MATLAB. Local asymptotic stability analysis was performed to investigate the models' behavior around equilibrium points. Moreover, the Lyapunov method was employed to establish global asymptotic stability of the models. Findings: The findings revealed insights into the dynamics of treating RAIR-DTC patients with tyrosine kinase inhibitors and immunotherapy. The stability analysis contributed to the understanding of the system's behavior under different conditions, offering valuable information for clinical application. Numerical simulations were then performed using Simbiology in MATLAB to simulate the treatment dynamics over time. Novelty: This research contributes to the field by presenting a novel mathematical model that integrates tyrosine kinase inhibitors and immunotherapy for the treatment of RAIR-DTC. The comprehensive analysis of equation stability enhances the reliability of the model, while the numerical simulations provide a practical and visual tool for assessing the potential outcomes of the proposed therapeutic strategy. This study offers a valuable framework for further exploration and refinement of treatment approaches for patients with challenging DTC cases. Keywords: Mathematical modeling, Routh Hurwitz Criteria, Jacobian, Differentiated thyroid Cancer, Stability Analysis