Estimating effects of drugs at different stages is directly proportional to the duration of recovery and duration of pulling through with the disease. For this reason, solving Pharmacokinetic models that investigate these effects is very important. In this study, numerical solutions of this type of one, two and three compartment nonlinear Pharmacokinetic models have been studied. Distributed order differential equations are used for the solution of the model. Numerical solutions have been found with the density function contained in distributed order differential equations and different values of this function. A Nonstandard finite difference scheme has been used for numerical solutions. Finally, stability analysis of equilibrium points of obtained discretized system has also been expressed with the help of the Schur-Cohn criteria.