To compute water-surface profiles in steady-state, gradually-varied flows in open-channel networks, an algorithm is presented which allows computation of the flow depths and discharges at all sections simultaneously instead of presently-used, step-by-step procedures. This algorithm is accurate, efficient, and suitable for computer solution, and it may be directly used for computing flows in parallel channels or in channel networks without necessitating graphical plots or other trial-and-error procedures. The Newton- Raphson method is used in the solution algorithm, and the resulting Jacobian matrix is transformed into a banded matrix in order to increase accuracy and to reduce computer time and storage. The governing equations are presented in detail and the solution algorithm is included in its entirety. Example problems are solved for illustration purposes and the results are compared with those obtained by the fourth-order Runge- Kutta method.