A numerical method of solution is developed for steady–state hyperbolic equations in conservation form. An upwind direction is identified and a smooth mapping is introduced to align the equations with that direction. A flux F(u) in the upwind direction is regarded as the ‘state’ variable and advanced numerically on a grid according to a second–order extension of Godunov' method. The solution u is determined from the computed flux, and a constraint on the mapping is that the Jacobian, det(∂F∂u), is of one sign. The method is applied to the equations modelling shock–wave propagation according to Whitham' theory of geometrical shock dynamics. In this context the solution vector u gives the components of the gradient of a function whose level curves determine the shock positions at successive times. Shock propagation in both uniform and non–uniform gases is considered and several example calculations are presented. A further application of the numerical method is inviscid supersonic flow, as governed by the steady Euler equations. The steady supersonic flow in a converging–diverging duct is computed to illustrate the numerical method for this latter application.