In this paper, 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 problem in the spectral domain to a quadratic optimization problem in the abundance map domain. In order to do so, Markov random field 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 and fully constrained least squares algorithms. Third, high resolution (HR) abundance maps are obtained using our proposed maximum a posteriori 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 data sets; namely the Cave, Harvard, and Hyperspectral Remote Sensing Scenes and compared with 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.