Recently, the simplified spherical harmonics equations (SPN) model has attracted much attention in modeling the light propagation in small tissue geometries at visible and near-infrared wavelengths. In this paper, we report an efficient numerical method for fluorescence molecular tomography (FMT) that combines the advantage of SPN model and adaptive hp finite element method (hp-FEM). For purposes of comparison, hp-FEM and h-FEM are, respectively applied to the reconstruction process with diffusion approximation and SPN model. Simulation experiments on a 3D digital mouse atlas and physical experiments on a phantom are designed to evaluate the reconstruction methods in terms of the location and the reconstructed fluorescent yield. The experimental results demonstrate that hp-FEM with SPN model, yield more accurate results than h-FEM with diffusion approximation model does. The phantom experiments show the potential and feasibility of the proposed approach in FMT applications.