  Morphoelasticity describes a family of nonlinear PDE system and has been
widely applied to understand growth and morphogenesis in living and nonliving
structures. Most previous studies depend on locally linearizing either the PDE
system or its discretized numerical counterpart. This involves conducting a
linear stability analysis for a trivial solution or utilizing Newton-Raphson
iteration on the force residual equations derived from finite element
approximation, respectively. The paper presents a numerical approach that
frames the search for solutions as the optimization of a discretized energy
functional, similar to Tallinen et al., PRL, 2013. In contrast to Tallinen et
al., we utilize the Hessian matrix of the energy to develop a gradient descent
algorithm that ensures energy stability, which allows us to search solutions
robustly in different types of bifurcations. The solutions are further refined
locally by Newton-Raphson iterations. We showcase the effectiveness of our
approach by studying shape transitions induced by the growth in the interior
layer of a bilayer annulus, with the outer layer constrained at its outermost
boundary. By investigating three cases with varying elastic-moduli ratios
between the two layers, we illustrate that our new method can not only identify
and characterize supercritical bifurcations leading to smooth wrinkles but also
capture subcritical transitions from smooth solutions to non-smooth
crease-forming solutions. Moreover, it reveals aperiodic solutions, both smooth
and non-smooth, that have not been previously reported.
