A novel amplitude factorization method is applied to solve a discrete buoyancy wave equation with arbitrary wind and temperature height distribution. The solution is given in the form of a cumulative product of complex factors, which are computed by a nonlinear, inhomogeneous, two-member recurrence formula, initiated from a radiative condition on top. Singularities of the wave equation due to evanescent winds are eliminated by turbulent friction. The method provides an estimation of the minimal vertical resolution, required to attain a stable accurate solution. The areas of application of the developed numerical scheme are high resolution modelling of orographic waves for arbitrary orography in general atmospheric stratification conditions, and testing of adiabatic kernels of numerical weather prediction models.