As of now, all analogue gravity models available in the literature deal with the emergence of an acoustic geometry through linear perturbations of transonic fluids only. It has never been investigated whether the analogue gravity phenomena is solely a consequence of linear perturbations, or rather a generic property of arbitrary perturbations of inhomogeneous, inviscid and irrotational fluids. In the present work, for the first time in the literature, we demonstrate that acoustic spacetimes may be formed through higher order non-linear perturbations, and thus establish that analogue gravity phenomena is rather more general than what was thought before. We consider spherically accreting astrophysical systems as a natural classical analogue gravity model, and develop a formalism to investigate non-linear perturbations of such accretion flows to arbitrary order. Our iterative approach involves a coupled set of equations for the mass accretion rate and the density of the fluid. In particular, we demonstrate that the wave equation for the mass accretion rate involves an acoustic metric which can be perturbatively constructed to all orders. We numerically solve the coupled equations about the leading transonic Bondi flow solution. This analysis uses boundary conditions set to the original unperturbed values, with the time dependence of the mass accretion rate perturbation taken to be exponentially damped. The perturbed solutions indicate that second order and higher perturbations of the metric generically cause the original acoustic horizon to oscillate and change in size. We explain this phenomenon in detail and its implications on non-linear perturbations of accretion flows in general.