The Bremmer series solution of the 3-D acoustical wave equation in generally inhomogeneous media, requires the introduction of pseudodifferential operators. These operators describe the (de)composition, propagation and interaction of the up- and downgoing waves. In this paper, a sparse matrix representation of the propagator is derived. In a similar way the matrix representations of the (de)composition and interaction of the up- and downgoing waves can be derived. The focus is on designing sparse matrices, keeping the accuracy high at the cost of ignoring any critical scattering-angle phenomena. Such matrix representations follow from well-known rational approximations of the vertical slowness, the Laplace operator and the vertical derivative. To reduce artificial azimuthal anisotropy a hexagonal prism-type grid is employed for the 3-D configuration. An optimization procedure is followed to minimize the errors in phase and group slowness, in the high-frequency limit for a given discretization rate. Through the Bremmer series it also accounts for the backscattered field. The resulting algorithm provides an improvement of the parabolic equation method. Numerical results will be presented at the conference.