In nonlinear conservation laws the flux function f(u) is usually single valued, but in many important applications it is hysteretic, i.e., it assigns different values depending on whether the input u(t) is increasing or decreasing in t. We present our recent results on a hysteresis model built with a collection of auxiliary ODEs under constraints. The model shares some similarities to, e.g., the well-known Preisach model, but is easy to calibrate and solve with either the resolvent or with the Yosida approximation to the constraint graphs. Due to the history dependence, the (stability) analysis is formulated on a product space rather than in (u) only. This is joint work with Ralph Showalter.