As an emerging and promising molecular imaging modality, bioluminescence tomography (BLT) can reconstruct the internal light source with the photon fluence on the small animal surface to reveal non-invasive molecular and cellular activities directly. In order to obtain higher precision and better spatial resolution in source reconstruction, the solution accuracy for the forward problem of BLT should be improved as high as possible. In this contribution, we present a modified element free Galerkin method (MEFGM) to calculate photon propagation in the biological tissue. This method is based on moving least squares (MLS) approximation which requires only a series of nodes in the region of interest, so complicated meshing task can be avoided compared with finite element method (FEM). Furthermore, MLS shape functions are further modified to satisfy the delta function property, which can simplify the processing of boundary conditions in comparison with traditional meshless methods. Finally, the numerical simulation experiments demonstrate the effectiveness and feasibility of this proposed method by comparing the solution of MEFGM with that of FEM.