We study the stability of the interface between (a) two adjacent viscous layers flowing due to gravity through an inclined or vertical channel that is confined between two parallel plane walls, and (b) two superimposed liquid films flowing down an inclined or vertical plane wall, in the limit of Stokes flow. In the case of channel flow, linear stability analysis predicts that, when the fluids are stably stratified, the flow is neutrally stable when the surface tension vanishes and the channel is vertical, and stable otherwise. This behaviour contrasts with that of the gravity-driven flow of two superimposed films flowing down an inclined plane, where an instability has been identified when the viscosity of the fluid next to the plane is less than that of the top fluid, even in the absence of fluid inertia. We investigate the nonlinear stages of the motion subject to finite-amplitude two-dimensional perturbations by numerical simulations based on boundary-integral methods. In both cases of channel and film flow, the mathematical formulation results in integral equations for the unknown interface and free-surface velocity. The properties of the integral equation for multi-film flow are investigated with reference to the feasibility of computing a solution by the method of successive substitutions, and a deflation strategy that allows an iterative procedure is developed. In the case of channel flow, the numerical simulations show that disturbances of sufficiently large amplitude may cause permanent deformation in which the interface folds or develops elongated fingers. The ratio of the viscosities and densities of the two fluids plays an important role in determining the morphology of the emerging interfacial patterns. Comparing the numerical results with the predictions of a model based on the lubrication approximation shows that the simplified approach can only describe a limited range of motions. In the case of film flow down an inclined plane, we develop a method for extracting the properties of the normal modes, including the ratio of the amplitudes of the free-surface and interfacial waves and their relative phase lag, from the results of a numerical simulation for small deformations. The numerical procedure employs an adaptation of Prony's method for fitting a signal described by a time series to a sum of complex exponentials; in the present case, the signal is identified with the cosine or sine Fourier coefficients of the interface and free-surface waves. Numerical simulations of the nonlinear motion confirm that the deformability of the free surface is necessary for the growth of small-amplitude perturbations, and show that the morphology of the interfacial patterns developing subject to finite-amplitude perturbations is qualitatively similar to that for channel flow.