In this study, we propose a novel single image Bayesian super-resolution (SR) algorithm where the hyperspectral image (HSI) is the only source of information. The main contribution of the proposed approach is to convert the ill-posed SR reconstruction (SRR) problem in the spectral domain to a quadratic optimization problem in the abundance map domain. In order to do so, Markov Random Field (MRF) based energy minimization approach is proposed and proved that the solution is quadratic. The proposed approach consists of five main steps. First, the number of endmembers in the scene is determined using virtual dimensionality. Second, the endmembers and their low resolution abundance maps are computed using simplex identification via the splitted augmented Lagrangian (SISAL) and fully constrained least squares (FCLS) algorithms. Third, high resolution (HR) abundance maps are obtained using our proposed maximum a posteriori (MAP) based energy function. This energy function is minimized subject to smoothness, unity and boundary constraints. Fourth, the HR abundance maps are further enhanced with texture preserving methods. Finally, HR HSI is reconstructed using the extracted endmembers and the enhanced abundance maps. The proposed method is tested on three real HSI datasets; namely the Cave, Harvard and Hyperspectral Remote Sensing Scenes (HRSS) and compared to state-of-the-art alternative methods using peak signal to noise ratio, structural similarity, spectral angle mapper and relative dimensionless global error in synthesis metrics. It is shown that the proposed method outperforms the state of the art methods in terms of quality while preserving the spectral consistency.