We investigate nonnegativity of exact and numerical solutions to a generalized Keller–Segel model. This model includes the so-called “minimal” Keller–Segel model, but can cover more general chemistry. We use maximum principles and invariant sets to prove that all components of the solution of the generalized model are nonnegative. We then derive numerical methods, using finite element techniques, for the generalized Keller–Segel model. Adapting the ideas in our proof of nonnegativity of exact solutions to the discrete setting, we are able to show nonnegativity of discrete solutions from the numerical methods under certain standard assumptions. One of the numerical methods is then applied to the minimal Keller–Segel model. Recalling known results on the qualitative behavior of this model, we are able to choose parameters that yield convergence to a nonhomogeneous stationary solution. While proceeding to exhibit these stationary patterns, we also demonstrate how naive choices of numerical methods can give physically unrealistic solutions, thereby justifying the need to study positivity preserving methods.