Summary: A robust method for solving the chemical non-equilibrium Navier-Stokes equations, including all of the species conservation and energy production equations, is developed. The algorithm is embodied in a fully coupled, implicit, large block structure. Van Leer flux splitting for inviscid terms and central differencing for viscous terms in the explicit operators are applied in the numerical algorithm. The fully coupled system is solved implicitly and the bi-conjugate gradient stable (Bi-CGSTAB) method with a preconditioner of incomplete lower-upper (LU)-factorization (ILU) is used for solving large block structure and diagonal dominate matrix equations. The computations are performed for the hypersonic flow over blunt bodies including half-cylinder, double ellipse and blunt nose. The adaptive grid constructed by moving grid method is employed to capture the shock location. Computational results in the present study are compared with other calculated data and exhibit good agreement. Convergence histories of the mean flow variables and species equations demonstrate that the fast convergence rate can be achieved by the preconditioned Bi-CGSTAB method.