  Entropy principles based on thermodynamic consistency requirements are widely
used for constitutive modeling in continuum mechanics, providing physical
constraints on a priori unknown constitutive functions. The well-known
M\"uller-Liu procedure is based on Liu's lemma for linear systems. While the
M\"uller-Liu algorithm works well for basic models with simple constitutive
dependencies, it cannot take into account nonlinear relationships that exist
between higher derivatives of the fields in the cases of more complex
constitutive dependencies.
  The current contribution presents a general solution set-based procedure,
which, for a model system of differential equations, respects the geometry of
the solution manifold, and yields a set of constraint equations on the unknown
constitutive functions, which are necessary and sufficient conditions for the
entropy production to stay nonnegative for any solution. Similarly to the
M\"uller-Liu procedure, the solution set approach is algorithmic, its output
being a set of constraint equations and a residual entropy inequality. The
solution set method is applicable to virtually any physical model, allows for
arbitrary initially postulated forms of the constitutive dependencies, and does
not use artificial constructs like Lagrange multipliers. A Maple implementation
makes the solution set method computationally straightforward and useful for
the constitutive modeling of complex systems.
  Several computational examples are considered, in particular, models of gas,
anisotropic fluid, and granular flow dynamics. The resulting constitutive
function forms are analyzed, and comparisons are provided. It is shown how the
solution set entropy principle can yield classification problems, leading to
several complementary sets of admissible constitutive functions; such problems
have not previously appeared in the constitutive modeling literature.
