  This article describes a numerical procedure designed to tune the parameters
of periodically-driven dynamical systems to a state in which they exhibit rich
dynamical behavior. This is achieved by maximizing the diversity of subharmonic
solutions available to the system within a range of the parameters that define
the driving. The procedure is applied to a problem of interest in computational
neuroscience: a circuit composed of two interacting populations of neurons
under external periodic forcing. Depending on the parameters that define the
circuit, such as the weights of the connections between the populations, the
response of the circuit to the driving can be strikingly rich and diverse. The
procedure is employed to find circuits that, when driven by external input,
exhibit multiple stable patterns of periodic activity organized in complex
tuning diagrams and signatures of low dimensional chaos.
