A two-stage source reconstruction algorithm for bioluminescence tomography (BLT) is developed using hybrid finite element method (FEM). The proposed algorithm takes full advantages of linear and quadratic FEMs, which can be used to localize and quantify bioluminescent source accurately. In the first stage, a large permissible region is roughly determined and then iteratively evolved to reduce matrix dimension using efficient linear FEM. In the final stage, high-convergence quadratic FEM is applied to improve reconstruction result. Both numerical simulation and physical experiment are performed to evaluate the proposed algorithm. The relevant results demonstrate that quantitative reconstruction can be well achieved in terms of computation efficiency, source position, power density, and total power when compared with previous studies. A two-stage source reconstruction algorithm for bioluminescence tomography (BLT) is developed using hybrid finite element method (FEM). The proposed algorithm takes full advantages of linear and quadratic FEMs, which can be used to localize and quantify bioluminescent source accurately. In the first stage, a large permissible region is significantly determined and then iteratively evolved to reduce matrix dimension using efficient linear FEM. In the final stage, high-convergence quadratic FEM is applied to improve reconstruction result. Both numerical simulation and physical experiment are performed to evaluate the proposed algorithm. The relevant results demonstrate that quantitative reconstruction can be well achieved in terms of computation efficiency, source position, power density, and total power when compared with previous studies.