In this work we introduce the finite volume (FV) approximation to the simplified spherical harmonics (SPN) equations for modeling light propagation in tissue. The SPN equations, with partly reflective boundary conditions, are discretized on unstructured grids. The resulting system of linear equations is solved with a Krylov subspace iterative method called the generalized minimal residual (GMRES) algorithm. The accuracy of the FV-SPN algorithm is validated through numerical simulations of light propagation in a numerical phantom with embedded inhomogeneities. We use a FV implementation of the equation of radiative transfer (ERT) as the benchmark algorithm. Solutions obtained using the FV-SPN (N > 1) algorithm are compared to solutions obtained with the ERT and the diffusion equation (SP1). Compared to the SP1, the SP3 solutions obtained using the FV-SPN algorithm can better approximate ERT solutions near boundary sources and in the vicinity of void-like regions. Solutions using the SP3 algorithm are obtained 9.95 times faster than solutions with the ERT-based algorithm.