A general method for computing solutions of systems of linear gravity-wave governing equations, allowing for vertically varying background parameters, is introduced and applied to two sets of governing equations for gravity waves in the thermosphere, both of which include molecular viscosity and thermal diffusion. The first set is fully compressible and the second set is anelastic. In the latter case, a new set of viscous and diffusive anelastic governing equations is given, and an anelastic dispersion relation, previously obtained only as a limiting case of a fully compressible dispersion relation, is shown to follow from the new anelastic governing equations. The new full-wave method is an extension of a numerical multilayer method described previously by Knight et al. (2019). Here, the term “numerical multilayer method” refers to a multilayer method that converges to a solution of an underlying vertical structure equation or system of equations in the limit of infinitesimal layer width. The previously described method was primarily suitable for dispersion-relation partial differential equations (PDEs), solutions of which only give approximate solutions of the underlying governing equations when the background parameters vary. Much of the theory developed for the previous method is still applicable to the new method, with some modifications, including the definition of upgoing and downgoing modes in terms of roots of the dispersion relation and the use of imaginary frequency shifts to make this classification of modes possible in all cases. The practical advantage of solving the anelastic equations instead of the compressible equations is that the needed imaginary frequency shifts are easier to determine in the anelastic case. To simplify the discussion, the application of the new full-wave method to the two sets of governing equations is only worked out for two-dimensional waves in the situation where kinematic viscosity varies but background temperature, scale height (for density), and horizontal wind are constant. Full-wave anelastic and compressible solutions are given for three idealized examples and compared with each other and with anelastic dispersion-relation PDE solutions. It is seen that full-wave anelastic solutions give good approximations of full-wave compressible solutions, while dispersion-relation PDE solutions become increasingly inaccurate with increasing kinematic viscosity for large horizontal wavelengths.