Internal gravity wave energy contributes significantly to the energy budget of the oceans, affecting mixing and the thermohaline circulation. Hence it is important to determine the internal wave energy flux J = p v, where p is the pressure perturbation field and v is the velocity perturbation field. However, the pressure perturbation field is not directly accessible in laboratory or field observations. Previously, a Green's function based method was developed to calculate the instantaneous energy flux field from a measured density perturbation field ρ(x, z, t), given a constant buoyancy frequency N . Here we present methods for computing the instantaneous energy flux J (x, z, t) for a spatially varying N (z), as in the oceans where N (z) typically decreases by two orders of magnitude from shallow water to the deep ocean. Analytic methods are presented for computing J (x, z, t) from a density perturbation field for N (z) varying linearly with z and for N 2 (z) varying as tanh(z). To generalize this approach to arbitrary N (z), we present a computational method for obtaining J (x, z, t). The results for J (x, z, t) for the different cases agree well with results from direct numerical simulations of the Navier-Stokes equations. Our computational method can be applied to any density perturbation data using the MATLAB graphical user interface EnergyFlux.