In this paper, we present a numerical study for HIV mathematical model of complex
order with medication resistance throughout therapy treatment. HIV is a virus that weakens
the immune system, making a person more susceptible to infections and diseases. This model consists of five nonlinear complex order differential equations where the derivatives is specified in the sense of Atangana-Baleanu-Caputo. Mittag-Leffler kernels are used in new numerical approaches to simulate complex order systems. These methods are based on Lagrange polynomial interpolation and the fundamental theorem of fractional calculus. For the two-step Lagrange polynomial interpolation, we suggest a straightforward adjustment to the step size to achieve high stability. The stability of the disease-free equilibrium point of the proposed model is presented. The complex order HIV model is mathematically studied using two different techniques: the standard and nonstandard Twostep Lagrange interpolation methods, which are suggested. To support the theoretical foundations, comparative investigations and numerical simulations are provided.