The K-distribution method is a fast approximative method used for calculating thermal infrared radiative transfer in the atmosphere, as opposed to the traditional Line-by-line method, which is precise, but very time-costly. Here we consider the atmosphere to consist of homogeneous and plane-parallel layers in local thermal equilibrium. This lets us use efficient upwards recursion for calculating the thermal infrared radiative transfer and ultimately the outgoing irradiance at the top of the atmosphere. Our specific implementation of the K-distribution method revolves around changing the integration space from the wavenumber domain to the g domain by employing Gauss-Legendre quadrature in two steps. The method is implemented in MATLAB and is shown to be several thousand times faster than the traditional Line-by-line method, with the relative error being only 3 % for the outgoing irradiance at the top of the atmosphere.