Numerical simulations of multi-scale flow problems such as hypersonic boundary layer transition, turbulent flows, computational aeroacoustics and other flow problems with complex physics require high-order methods with high spectral resolutions. For instance, the receptivity mechanisms in the hypersonic boundary layer are the resonant interactions between forcing waves and boundary-layer waves, and the complex wave interactions are difficult to be accurately predicted by conventional low-order numerical methods. High-order methods, which are robust and accurate in resolving a wide range of time and length scales, are required. The objective of this dissertation is to develop and analyze new very high-order numerical methods with spectral-like resolution for flow simulations on structured grids, with focus on smooth flow problems involving multiple scales. These numerical methods include: the multi-layer compact (MLC) scheme, the directional multi-layer compact (DMLC) scheme, and the least square multi-layer compact (LSMLC) scheme.In the first place, a new upwind multi-layer compact (MLC) scheme up to seventh order is derived in a finite difference framework. By using the ‘multi-layer’ idea, which introduces first derivatives into the MLC schemes and approximates the second derivatives, the resolution of the MLC schemes can be significantly improved within a compact grid stencil. The auxiliary equations are introduced, and they are the only nontrivial equations. The original equation requires no approximation which contributes to good computational efficiency. In addition, the upwind MLC schemes are derived on centered stencils with adjustable parameters to control the dissipation. Fourier analysis is performed to show that the new MLC schemes have very small dissipation and dispersion in a wide range of wavenumbers in both one and two-dimensional cases, and the anisotropic error is much smaller than conventional finite difference methods in the two-dimensional case. Comparison with discontinuous-Galerkin methods is performed with Fourier analysis as well. Furthermore, stability analysis with matrix method shows that high-order boundary closure schemes are stable because of compactness of the stencils. The accuracies and rates of convergence of the new schemes are validated by numerical experiments of the linear advection equation, the nonlinear Euler equations, and the Navier-Stokes equations in both one and two-dimensional settings. The numerical results show that good computational efficiency, very high-order accuracies, and high spectral resolutions especially on coarse meshes can be attained with the MLC scheme. On the other hand, even though the MLC scheme is promising in most test cases, it shows weak numerical instabilities for a small range of wavenumbers when it is applied to multi-dimensional flows, which are mainly triggered by the inconsistency between its one and two-dimensional formulations. The instability could lead to divergence in long-time multi-dimensional simulations. Moreover, the cross-derivative approximation in the MLC scheme requires an ad-hoc selection of supporting grid points, and the cross-derivative approximation is relatively inefficient for very high-order cases. To address the remaining challenges of the MLC scheme and achieve better performance for multi-dimensional flow simulations, another two new schemes are developed – the directional multi-layer compact (DMLC) scheme, and the least square multi-layer compact (LSMLC) scheme. In the second place, a new upwind directional multi-layer compact (DMLC) scheme is developed for multi-dimensional simulations. The main idea of the DMLC scheme is to introduce auxiliary equation for cross derivative in multi-dimensional cases. Consequently, the spatial discretization can be fulfilled along each dimension independently. With this directional discretization technique, the one-dimensional formulation of the MLC scheme can be applied to all spatial derivatives in a multi-dimensional governing equation. Therefore, the DMLC scheme overcomes the inconsistency between one and two-dimensional formulations of the MLC scheme, and it also avoids the ad-hoc cross-derivative approximations. Two-dimensional Fourier analysis demonstrates that all modes of the DMLC scheme are stable in the full range of wavenumbers, and it has better spectral resolution and smaller anisotropic error than the MLC scheme. Stability analysis with matrix method indicates that stable boundary closure schemes are much easier to be obtained in the DMLC scheme. Numerical tests in the linear advection equation and the nonlinear Euler equations validate that the DMLC scheme are more accurate and require less CPU time than the MLC scheme on the same mesh. In particular, the long-time simulation results reveal that the DMLC scheme is always stable for both periodic and non-periodic boundary conditions in two-dimensional cases. In the third place, a new upwind least square multi-layer compact (LSMLC) scheme is developed for multi-dimensional simulations. The main idea of the LSMLC scheme is using the weighted least square approximation to redesign the two-dimensional formulation for cross derivatives. It avoids the ad-hoc selection of grid points in the MLC scheme. Meanwhile, the two-dimensional upwind scheme can be derived by introducing upwind correction into the weight function. The upwind factor β can adjust the dissipation and stability of the LSMLC scheme. Lagrange multiplier is used to ensure that the LSMLC scheme satisfies both the consistency constraint at the base point and the one-dimensional constraint from the MLC scheme. The LSMLC scheme does not increase computational cost on structured meshes, and can be implemented in the same way as the MLC scheme. A parametric study based on two-dimensional Fourier analysis shows that the truncated Gaussian distribution (TGD) weight function leads to better LSMLC scheme among other weight functions because it removes the numerical instability and maintain small dissipations. The LSMLC scheme has larger dissipation than the MLC scheme, and shows similar spectral resolution. Stability analysis with matrix method indicates that a combination of an interior LSMLC scheme and MLC boundary closure schemes can improve the boundary stability while maintaining small dissipation. Numerical tests in the linear advection equation and the nonlinear Euler equations validate that the LSMLC scheme produces slightly larger errors compared with MLC scheme. The long-time simulation results reveal that the LSMLC scheme is always stable for both periodic and non-periodic boundary conditions in two-dimensional cases.Overall, the new very high-order multi-layer compact finite difference methods have the properties of simple formulations, high-order accuracies, spectral-like resolutions, and compact stencils, and they are suitable for accurate simulation of smooth multi-scale flows with complex physics. Among the three schemes developed in this dissertation, the DMLC scheme is always the best choice for multi-dimensional simulations because it shows comprehensive improvements from the MLC scheme with consistent stability, higher accuracy and spectral resolution, and better computational efficiency. The LSMLC scheme is also appropriate considering it has consistent stability and it is easy to be implemented.