In the lectures, I will first introduce and review several chemotaxis models including the classical Patlak-Keller-Segel model. Chemotaxis is the phenomenon in which cells, bacteria, and other single-cell or multicellular organisms direct their movements according to certain chemicals (chemoattractants) in their environment. Chemotaxis is an important process in many medical and biological applications including bacteria/cell aggregation, pattern formation mechanisms, and tumor growth. The mathematical models of chemotaxis are usually described by highly nonlinear time dependent systems of partial differential equations (PDEs). Therefore, accurate and efficient numerical methods are very important for the validation and analysis of these systems. Furthermore, a common property of all existing chemotaxis systems is their ability to model a concentration phenomenon that mathematically results in solutions rapidly growing in small neighborhoods of concentration points/curves. The solutions may blow up or may exhibit a very singular, spiky behavior. This blow-up represents a mathematical description of a cell concentration phenomenon that occurs in real biological systems. In either case, capturing such solutions numerically is a challenging problem. In the lectures, we will conduct a detailed review of several recently developed numerical methods for the approximation and simulation of the chemotaxis and related models. In particular, we will discuss several finite element, finite-volume and hybrid finite-volume-finite-difference methods, which are designed to achieve an ultimate goal of developing a highly accurate, stable and robust numerical methods for chemotaxis models.