A fractional order mathematical model that already exists in the literature, was considered. This model was established to study the effects of medicinal treatment in people infected with the human immunodeficiency virus (HIV). The importance of this study is that the model evaluates, among other parameters, the density of healthy and HIV-infected CD4+ T cells. These data are very necessary for the subject infected by the virus given the effects that an antiretroviral treatment causes in it. The objective of this work is to consider several numerical schemes that involve fractional derivatives in order to compare their behaviors and to obtain a good approximation of the mentioned model solution. Convergence of these schemes will be studied as well as sensitivity with respect to the variation of the parameters η (drug efficacy) and α (fractional derivative order). Furthermore, through the collection of medical records of people living with HIV, it is intended to determine the optimal fractional derivative order for the model and to compare it with the classical model.