This thesis is devoted to developing robust and efficient high-order finite volume solvers for numerically solving the two-dimensional steady Euler equations. The associated numerical framework consists of a Newton iteration method to linearize the non-linear governing equations, and a geometrical multigrid method to solve the derived linear system, In the simulation, instead of using the standard regularization where a time-stepping relaxation term is added into the Jacobian matrix of the Newton-iteration, the l¹ norm of the local residual of each element is used to regularize the Jacobian matrix of the Newton iteration. In order to obtain finite volume schemes with second and third order of accuracy on unstructured grids, the integrated linear reconstruction(ILR)and the non-oscillatory k-exact reconstruction are used for the solution reconstruction, respectively, The IR was originally proposed for the time-dependent conservation laws, and the numerical results show that the system residual cannot converge smoothly to steady state with the direct use of ILR. To resolve the issue, two simple gradient smoothing procedures are proposed. The convergence to steady state by using the smoothed gradients is verified by various numerical examples. On the other hand, to obtain physically correct numerical solutions, the Non-Uniform Rational B-Splines (NURBS) are used to describe the curved boundary. The high-order numerical methods have been widely studied in the research field of computational fluid dynamics (CFD). In this thesis, the non-oscillatory k-exact reconstruction is used to obtain high-order finite volume scheme. To obtain a high quality high-order numerical scheme and to obtain physically correct numerical results, the NURBS are also used to describe the curved boundary. The accuracy as well as the convergence rate of the numerical methods can be improved by using NURBS to describe the curved boundary. Several numerical examples are presented to demonstrate the effectiveness of using the non-oscillatory k-exact reconstruction with NURBS. To further improve the efficiency of the high-order finite volume scheme, the following two strategies are considered in this thesis. First, a new reconstruction patch is specially designed for the boundary elements, which could significantly improve the convergence to the steady state by using the Newton-type fnite volumemethod. Second, the adjoint-based a posteriori error estimate is used to design theerror indicator for the h-adaptive method with the aim of reducing the error in thequantity of interest. Finally, the conclusion of this thesis and some possible future research topics will be given in the last part.