A mathematical model of the hypothalamic-pituitary-adrenal (HPA) axis of the human endocrine system is proposed. This new model provides an improvement over previous models by introducing two nonlinear factors with physiological relevance: 1) a limit to gland size; 2) rejection of negative hormone concentrations. The result is that the new model is by far the most robust; e.g., it can tolerate at least -50% and +100% perturbations to any of its parameters. This high degree of robustness allows one, for the first time, to model features of the system such as circadian rhythm and response to hormone injections. In addition, relative to its closest predecessor, the model is simpler; it contains only about half of the parameters, and yet achieves more functions. The new model provides opportunities for teaching endocrinology within a biological or medical school context; it may also have applications in modeling and studying HPA axis disorders, for example, related to gland size dynamics, abnormal hormone levels, or stress influences.