Includes bibliographical references.Future large-scale structure surveys are expected to pin-down the properties of dark energy significantly more by mapping the cosmic web to unprecedented precision. To take advantage of such state-of-the-art technologies, the evermore accurate modelling of structure formation is absolutely necessary. While relativistic linear and non-relativistic (Newtonian) non-linear effects have been well established (although improvements are still being made), a fairly unexplored area is the impact of relativistic, non-linear effects on structure formation. As an attempt in this direction, we consider linear perturbations of a Lemaître-Tolman-Bondi (LTB) spacetime. LTB models are spherically symmetric but inhomogeneous exact dust solutions to the Einstein field equations. They are known to accommodate most observations of the background universe without dark energy. In this work we present a new numerical code to solve the set of coupled partial differential equations that describe the evolution of the (polar) perturbations, test it in the case of a Hubble-scale LTB void, and demonstrate its excellent stability and convergence. We then explore the solutions for a variety of generic initial conditions. The variable that closely resembles the Newtonian potential is shown to excite propagating (tensor) as well as rotational (vector) modes at the percent-level. Comparing our results to that which ignores the full coupling, we estimate percent-level corrections to the amplitude of the galaxy correlation function when only the scalar degrees of freedom are included. In addition, we showed that the anisotropic correlation function can nevertheless be used as a test of the Copernican Principle. Note that our code has applications to other scenarios as well in which spherical symmetry is a good approximation, such as the lensing of gravitational waves by intervening halos/voids
