A new method of characteristics numerical scheme for three- dimensional steady flow was developed which has second-order accuracy. Heretofore all such schemes for three-dimensional flow have had accuracies less than second-order. A complete numerical algorithm for computing internal supersonic flows of the type encountered in ramjet, scramjet or rocket propulsion systems was developed and programmed for both the IBM 7094 and CDC 6500 computers. The method was tested for order of accuracy using the exact solution for source flow and Prandtl-Meyer flow. The results of these tests verified the second-order accuracy of the scheme. Additional accuracy tests using existing methods for solution of two-dimensional axisymmetric flows showed that the scheme produces accuracies comparable to that of the two-dimensional method of characteristics. The computer program was used to generate the flow field for several three-dimensional nozzle contours and for nonsymmetric flow into an axisymmetric nozzle. These results reveal the complex nature of three- dimensional flows and the general inadequacy of quasi-three-dimensional analyses which neglect crossflow. An operationally convenient computer program was produced. The program has the capability to analyze nonisoenergetic and nonhomentropic flows of a calorically perfect gas or homentropic flows of a real gas in chemical equilibrium. The initial-value surface options include uniform flow, source flow or axisymmetric tabular data. The nozzle boundary options include conical nozzles, axisymmetric contoured nozzles and super-elliptical nozzles.