We use one-dimensional numerical simulations to study spherical collapse in the ƒ(R) gravity models. We include the nonlinear self-coupling of the scalar field in the theory and use a relaxation scheme to follow the collapse. We find an unusual enhancement in density near the virial radius which may provide observable tests of gravity. We also use the estimated collapse time to calculate the critical overdensity δc used in calculating the mass function and bias of halos. We find that analytical approximations previously used in the literature do not capture the complexity of nonlinear spherical collapse.