The strain image contrast of some in vivo breast lesions changes with increasing applied load. This change is attributed to differences in the nonlinear elastic properties of the constituent tissues suggesting some potential to help classify breast diseases by their nonlinear elastic properties. A phantom with inclusions and long-term stability is desired to serve as a test bed for nonlinear elasticity imaging method development, testing, etc. This study reports a phantom designed to investigate nonlinear elastic properties with ultrasound elastographic techniques. The phantom contains four spherical inclusions and was manufactured from a mixture of gelatin, agar and oil. The phantom background and each of the inclusions have distinct Young's modulus and nonlinear mechanical behavior. This phantom was subjected to large deformations (up to 20%) while scanning with ultrasound, and changes in strain image contrast and contrast-to-noise ratio between inclusion and background, as a function of applied deformation, were investigated. The changes in contrast over a large deformation range predicted by the finite element analysis (FEA) were consistent with those experimentally observed. Therefore, the paper reports a procedure for making phantoms with predictable nonlinear behavior, based on independent measurements of the constituent materials, and shows that the resulting strain images (e.g., strain contrast) agree with that predicted with nonlinear FEA.