Newton iteration, in combination with the conjugate gradient algorithm Generalized Minimal Residual (GMRES), forms a powerful tool for implicit numerical time integration of nonlinear partial differential equations such as the time-dependent Navier–Stokes equations. Particularly for steady compressible fluid flow, it allows large time steps to faster attain the asymptotical steady state of flow on a computer. Due to the necessity to resolve all relevant flow phenomena, the systems of linear equations involved may thereby become very large, so that the storage requirements of the matrices may surmount the available computer memory. There are two other disadvantages if nonlinearity prevents complete linearization. Then firstly, the desired quadratic convergence behaviour of Newton’s method is lost, and secondly, divergence can occur if the flow problem contains strong discontinuities like shocks. Looking at GMRES, there is no operation involving the matrix explicitly. Actually, only inner products with vectors must be formed, which represent weighted derivatives in certain directions and can be approximated consequently in the sense of finite differences. This way, no linearization remains necessary, and the memory requirement of the code decreases dramatically, in particular for multi-dimensional problems. Based on these arguments, the paper presents a finite difference version of GMRES (FD-GMRES) which considerably improves the convergence and reduces the required computer memory. The numerical errors caused by finite precision arithmetics and by truncation of the Taylor series are evaluated in competition to get an optimized error bound for FD-GMRES. With the help of test cases such as the flow across a cylinder (2-d) in the hypersonic regime, we show that this method is superior to the standard algorithm using the Fréchet Differentials explicitly. The improved stability of the method results in much faster convergence.