We demonstrate the modeling of optical response of exciton-polarons based on the well established Holstein Hamiltonian to model coupled exciton-phonon systems in organic molecular chains. Our approach uses Green’s functions to compute the density of states and the linear optical susceptibility, and thus eliminates the conventional and computationally expensive step of diagonalizing a large Hamiltonian matrix. We exploit this technique further to focus exclusively on the optically active states when computing the linear optical response, and significantly reduce the computational effort to construct the optical susceptibility. In this article, we demonstrate the computation of absorption and emission spectra of Alq3 at 4.2 K and at room temperature using our model. Using the two parameters of the Holstein model, the inhomogeneous broadening energies, and a phenomenological reorganization energy of the solute, we obtain excellent fits to established experimental results. We then use this model inside the larger simulation of a 3-layer organic light emitting (OLED) structure composed of Alq3, Alq3:DCJTB, and α-NPD, which are the electron transport, emissive, and hole transport layers respectively. In our methodology, we also couple the optical response into the rate equations for exciton dynamics in addition to computing the spectrum of light output by the device.