A mixed finite element method combining an iso-parametric $Q_2$-$P_1$ element and an iso-parametric $P_2^+$-$P_1$ element is developed for the computation of multiple cavities in incompressible nonlinear elasticity. The method is analytically proved to be locking-free and convergent, and it is also shown to be numerically accurate and efficient by numerical experiments. Furthermore, the newly developed accurate method enables us to find an interesting new bifurcation phenomenon in multi-cavity growth.