The human brain is characterized by its folded structure, being the most folded brain among all primates. The process by which these folds emerge, called gyrogenesis, is still not fully understood. The brain is divided into an outer region, called gray matter, which grows at a faster rate than the inner region, called white matter. It is hypothesized that this imbalance in growth -- and the mechanical stress thereby generated -- drives gyrogenesis, which is the focus of this thesis. Finite element simulations are performed where the brain is modeled as a non-linear elastic and growth is introduced via a multiplicative decomposition. A small section of the brain, represented by a rectangular slab, is analyzed. This slab is divided into a thin hard upper layer mimicking the gray matter, and a soft substrate, mimicking the white matter. The top layer is then grown tangentially, while the underlying substrate does not grow. JuFold, the software developed to perform these simulations, is introduced, and its design is explained. An overview of its capabilities, and examples of simulation possibilities are shown. Additionally, one patent-leading application of JuFold in the realm of material science showcases its flexibility. Simulations are first performed by minimizing the elastic energy, corresponding to the slow growth regime. Systems with homogeneous cortices are studied, where growth initially compresses, and then buckles the cortical region, which generates wavy patterns with wavelength proportional to cortical thickness. After buckling, the sulcal regions (i.e. the valleys of the system) are thinner than the gyral regions (i.e. the hills). Introducing thickness inhomogeneities along the cortex lead to new and localized configurations, which are strongly dependent not only on the thickness of the region, but also on its gradient. Furthermore, cortical landmarks appear sequentially, consistent with the hierarchical folding observed during gestation. A linear stability theory is developed based on thin plate theory and is compared with homogeneous and inhomogeneous systems. Next, we turn to more physically stringent dynamic simulations. For slow growth rate and time-constant thickness, the results obtained through energy minimization are recovered, justifying previous literature. For faster growth, an overshoot of the wavenumber and a broad wavenumber spectrum are observed immediately after buckling. After a relaxation period, where the average wavenumber decreases and the wavenumber spectrum narrows, it is observed that the system stabilizes into a finite spectrum, whose average wavelength is smaller than that expected from energy minimization arguments. Cortical inhomogeneities are further explored in this new regime. Systems with inhomogeneous cortical thickness are revisited, with effects similar to the homogeneous cortex (i.e., results are consistent between the slow growth and the quasistatic regimes, and overshoot is observed in the fast growth regimes). Systems with inhomogeneous cortical growth are simulated, with this new type of inhomogeneity inducing fissuration and localized folding. The interplay between these two inhomogeneities is studied, and their interaction is shown to be nonlinear, with each inhomogeneity type inhibiting the folding effects of the other. That is, the folding profile of each individual region emerges as a result of the local inhomogeneity, and the system does not display an intermediate behavior. Finally, these results are compared with an extended linear stability theory. Taken together, our simulations and analytical theory expose new phenomena predicted by an incremented buckling hypothesis for folding and show a series of new avenues which could give rise to the important cortical features in the mammalian brain, especially those related to higher-order folding.