Growth-elasticity (also known as morphoelasticity) is a powerful model framework for understanding complex shape development in soft biological tissues. At each instant, by mapping how continuum building blocks have grown geometrically and how they respond elastically to the push-and-pull from their neighbors, the shape of the growing structure is determined from a state of mechanical equilibrium. As mechanical loads continue to be added to the system through growth, many interesting shapes, such as smooth wavy wrinkles, sharp creases, and deep folds, can form on the tissue surface from a relatively flatter geometry.Previous numerical simulations of growth-elasticity have reproduced many interesting shapes resembling those observed in reality, such as the foldings on mammalian brains and guts. In the case of mammalian guts, it has been shown that wavy wrinkles, deep folds, and sharp creases on the interior organ surface can be simulated even under a simple assumption of isotropic uniform growth in the interior layer of the organ. Interestingly, the simulated patterns are all regular along the tube’s circumference, with either all smooth or all sharp indentations, whereas some undulation patterns in reality exhibit irregular patterns and a mixture of sharp creases and smooth indentations along the circumference. Can we simulate irregular indentation patterns without further complicating the growth patterning?In this paper, we have discovered abundant shape solutions with irregular indentation patterns by developing a Rayleigh–Ritz finite-element method (FEM). In contrast to previous Galerkin FEMs, which solve the weak formulation of the mechanical-equilibrium equations, the new method formulates an optimization problem for the discretized energy functional, whose critical points are equivalent to solutions obtained by solving the mechanical-equilibrium equations. This new method is more robust than previous methods. Specifically, it does not require the initial guess to be near a solution to achieve convergence, and it allows control over the direction of numerical iterates across the energy landscape. This approach enables the capture of more solutions that cannot be easily reached by previous methods. In addition to the previously found regular smooth and non-smooth configurations, we have identified a new transitional irregular smooth shape, new shapes with a mixture of smooth and non-smooth surface indentations, and a variety of irregular patterns with different numbers of creases. Our numerical results demonstrate that growth-elasticity modeling can match more shape patterns observed in reality than previously thought.