In the search for causes and cures of cancer diseases, many mathematical models developed have resulted in systems of nonlinear stiff ordinary differential equations. With these models, many numerical estimates of biological knowledge of the parameters have been obtained, a number of phenomena interpreted, and predictions were made in order to gain further knowledge of cancer development and possible treatment. In this study, numerical simulations of the models were performed using continuous block implicit hybrid methods and the results obtained support the theoretical and clinical findings. We analyzed the interactions among the various tumor cell populations and present the results graphically. From the graphical representation of results, one can clearly see the effects of all the tumor cell populations involved in the competition, as well as the effects of some treatments by the applications of some therapeutic agents which have been heavily used in the clinical treatments of breast cancer. The treatments in the past were mostly conventional chemotherapies, which were used either singly (alone) or in combination with each other or other therapies, and all played vital roles, except for the side effects that these therapies incur in normal tissues and organs. Thus, from recent research works, it is now clear that in many cases they do not represent a complete cure. Therefore, the need to address not only the preventative measures of breast cancer, but also more successful treatment, is clear, and can be successfully achieved to increase the survival rate of breast cancer patients.