<p>The emergence of drug resistance remains a major challenge in cancer treatment, necessitating the development of mathematical tools to explore effective therapeutic strategies. In this paper, we present a novel mathematical model governed by a bilinear system of ordinary differential equations, designed to capture the dynamic interactions between two groups of cancer subpopulations under combined therapies. This model is distinctive in that it integrates two key approaches in cancer modeling drug resistance dynamics and cell cycle-specific progression allowing us to study both how cells acquire resistance and how they progress through different stages of growth. The model incorporates phenotypic transitions, proliferation across cell cycle phases, and therapeutic responses. Through a detailed stability analysis, we derive critical threshold values that determine the persistence or elimination of the tumor population, and rigorously establish the existence and nature of the system’s equilibria. To investigate treatment optimization, we formulate and analyze an optimal control problem aimed at minimizing the resistant cell population while limiting treatment intensity and preserving the sensitive cell line. Using Pontryagin’s Maximum Principle, we derive the necessary conditions for optimality and employ a numerical forward–backward sweep method to solve the resulting system. Numerical simulations implemented in MATLAB illustrate how specific treatment combinations and scheduling can significantly influence tumor progression and therapeutic outcomes. The results highlight the potential of mathematical modeling to guide the design of combination therapies and support clinical decision-making.</p>