We explore the use of PDE constrained nonlinear optimization techniques to optimize and design electromagnetic crystals which exhibit frozen mode behavior. This is characterized by Van Hove singularities in the dispersion relation, e.g., stationary reflection points and degenerate band edge points. Hence, the optimization process modifies the dispersion relation by adjusting the geometries and material parameters. The resulting algorithm is found to be capable of recovering all known crystal configurations as well as many new configurations, some of which display dramatically improved properties over previously used configuration. We investigate both gyrotropic photonic crystals and degenerate band edge crystals as well as the more complex case of the oblique incidence. In this latter case, we extend the investigation to the three-dimensional case to identify the first three-dimensional crystal exhibiting frozen mode behavior.