The objective of this paper is to provide a unified treatment of elastic and electromagnetic (EM) wave propagation in horizontally layered media for which the parameters in the partial differential equations are piece-wise continuous functions of only one spatial variable. By applying a combination of Fourier, Laplace, and Bessel transforms to the partial differential equations describing the elastic or EM wave propagation I obtain a system of 2n linear ordinary differential equations. The 2n X 2n coefficient matrix is partitioned into 4 n X n submatrices. By a proper choice of variables, the diagonal submatrices are zero and the off-diagonal submatrices are symmetric.All the results in the paper are derived from the symmetry properties of this general equation. In the appendices it is shown that three-dimensional elastic waves, cylindrical P-SV waves, acoustic waves, and electromagnetic waves in isotropic layered media can all be represented by an equation with the same properties.The symmetry properties of the system matrix are used to derive simplified equations for computing the propagator matrix for a stack of inhomogeneous layers. The wave field is also decomposed into upgoing and downgoing waves by an eigenvector decomposition which is much simplified compared with the general case of a full 2n X 2n system matrix. This wave field decomposition is a suitable starting point for deriving one-way wave equations and WKB-approximations of different order. These approximations have potential applications in generalized migration schemes.Two propagation invariants are derived using the symmetry properties of the system matrix. One of these is only valid for lossless media and corresponds to the conservation of energy.For a stack of inhomogeneous layers, transmission and reflection matrices for upward and downward propagation are defined. Using the two propagation invariants, a number of symmetry properties of the reflection and transmission matrices are derived. The relationship between the reflection and transmission matrices and the propagator matrix is given. Computation of the reflection and transmission matrices for two inhomogeneous layers is done by Redheffer's star product. This composition rule has been derived for P-SV waves by Kennett. The inverse of the star product, apparently unknown in seismology, is also given. This is a rule which may be used to remove the effect of an inhomogeneous layer at the top or bottom of a stack of layers. Such layer stripping techniques have possible applications in general inversion schemes. It is also shown that the reflection and transmission matrices of an inhomogeneous medium can be found by solving a matrix Riccati equation.For a stack of inhomogeneous layers bounded above by a free surface, modified reflection and transmission matrices are defined. Using the two propagation invariants, a number of symmetry properties of the modified reflection and transmission matrices are derived. For lossless media a generalized Kunetz equation is given. The modified reflection and transmission matrices are expressed in terms of the partitioned submatrices of the propagator matrix and in terms of the usual reflection and transmission matrices.I also derive the response of a buried point source for a layered medium bounded by a free surface and a homogeneous half-space, and for a layered medium bounded by two homogeneous half-spaces. More general sources can be treated by superposition. In an appendix, I use multidimensional Fourier transforms to derive the decomposition of a spherical wave into plane waves (the Weyl integral) and into cylindrical waves (the Sommerfeld integral). Both these decompositions include inhomogeneous or evanescent waves. The Whittaker integral represents a decomposition into traveling waves only, but an explicit form of the Whittaker integral for spherical waves does not appear to be known.The decomposition into upgoing and downgoing waves breaks down for horizontally traveling waves such as channel waves and surface waves. The reflection and transmission matrices do not exist in this case. This fact is used to derive dispersion relationships for channel waves and surface waves by requiring that certain determinants shall be zero.