The deviations of non-linear perturbations of black holes from the linear case are important in the context of ringdown signals with large signal-to-noise ratio. To facilitate a comparison between the two, we derive several results of linear perturbation theory in coordinates which may be adopted in numerical work. Firstly, we address the questions: for an initial configuration of a massless scalar field, what is the amplitude of the excited quasinormal mode (QNM) for any observer outside outside the event horizon, and furthermore what is the resulting tail contribution? This is done by constructing the full Green’s function for the problem with exact solutions of the confluent Heun equation satisfying appropriate boundary conditions. We then present an implementation of the dual foliation generalized harmonic gauge (DF-GHG) formulation within our pseudospectral code 'bamps'. The formalism promises to give greater freedom in the choice of coordinates that can be used in numerical relativity. As a specific application we focus here on the treatment of black holes in spherical symmetry. Existing approaches to black hole excision in numerical relativity are susceptible to failure if the boundary fails to remain outflow. We present a method, called DF-excision, to avoid this failure. We compare the results of DF-excision with a naive setup and show that DF-excision proves reliable even when the previous approach fails. Finally, we use DF-excision and DF-GHG to perform fully non-linear simulations of a massless scalar field minimally coupled to general relativity using coordinates which are analogous to the Kerr-Schild coordinates from the linear setting. We perform a detailed analysis comparing the results of the non-linear simulations and the linear simulations and perform modeling of the non-linear time series data using several machine learning algorithms and compare their performance against each other.