High order numerical methods present the desirable property of having much less error than second order methods for the same number of nodes, or degrees of freedom. This makes them very attractive for applications like computational acoustics and aeroacoustics, prediction of transition or direct numerical simulation of turbulence. But the fact that high order numerical methods are prone to exhibit numerical instabilities have hindered their widespread use. These instabilities have been related to the presence of boundaries and the way in which boundary conditions are implemented in these methods. In the present work numerical stability of high order finite difference methods is achieved by suppression of the Runge phenomenon responsible for the vigorous oscillations of polynomial interpolations performed on equispaced nodes. By following the philosophy behind the Chebyshev polynomials, the idea of an optimum grid for piecewise polynomial interpolations of degree q _ N is introduced. It is then shown, that when q = N, this polynomial interpolation on optimum nodes coincides with the Chebyshev interpolation, and that the resulting finite difference schemes are therefore equivalent to Chebyshev collocation methods. In the opposite limiting case, when q _ N, it is shown that only a few nodes, of O(q), close to the boundaries need to be clustered to achieve the desired stability of the numerical method. Through numerical solutions of the convection-diffusion operator and the wave operator it is shown, that this new family of finite difference methods achieves numerically stable solutions for any degree q < N and that these solutions present correct transient behaviours. This has motivated the use of these numerical methods for the study of two different fluid-dynamic combustion problems. In the first problem the local flame extinction and re-ignition of a counterflow diffusion flame perturbed by a laminar vortex ring is investigated. The local extinction of the flame leads to the apparition of flame-edges separating the burning and extinguished regions of the perturbed mixing layer. The dynamics of these edges is modelled using previous numerical results, where heat release effects are fully taken into account, that provide the propagation velocity of triple- and edge-flames in terms of the local Damk¨ohler number. The temporal evolution of the mixing layer is determined using the classical mixture fraction approach, assuming a simplified analytical description of the flow field, where both unsteady and curvature effects are taken into account. Although variable density effects play an important role in exothermic reacting mixing layers, in this work the description of the mixing layer is carried out using the constant density approximation. The theoretical model reveals the relevant non-dimensional parameters governing diffusion-flame/vortex interactions, and provides the parameter range for the more relevant regime of local flame extinction followed by re-ignition via flame-edges. Despite the simplicity of the model, the results show very good agreement with previously published experimental results of diffusion-flame/vortex interactions. The second fluid-dynamic combustion problem deals with the vaporization and subsequent combustion of fuel droplets immersed in slowly convective flows. For small values of the Peclet number Pe, the convection associated with the velocity of the oxidizer stream is in first approximation negligible at distances to the droplet of the order of the droplet radius a. Only in the Oseen region, located at distances of the order of a/Pe, do the convective effects become as important as the diffusive ones. For typical hydrocarbon fuels, where the overall stoichiometric ratio S is large compared to unity, the flame is located in the very same region if the distinguished limit of Pe _ 1/S is considered, inducing temperature and density variations of order unity, which require the use of numerical techniques for the description of the resulting fluid-dynamic problem. The overall analysis of this multiscale problem is carried out using matched asymptotic expansions between the different regions of the flow, where the matching has to be done between the semi-analytical solutions obtained for the inner region and the numerical solutions obtained for the Oseen region. The presented analysis reveals the non-dimensional parameters that are relevant in each of the regions and shows that the presence of the flame significantly modifies the convective velocity felt by the droplet, thus altering its vaporization rate and drag.