This talk presents a numerical analysis for quantitative photoacoustic tomography, aiming to reconstruct optical coefficients (diffusion and absorption) using internal data. Our approach solves an inverse diffusivity problem and an elliptic direct problem. The stability of the inverse problem significantly depends on a non-zero condition in the internal observations, a condition that can be met using randomly chosen boundary excitation data. Utilizing these randomly generated boundary data, we implement an output least squares formulation combined with finite element discretization to solve the inverse problem. We provide a rigorous error estimate in $L^2(\\Omega)$ norm for the numerical reconstruction using a weighted energy estimate, motivated by the analysis of a newly proposed conditional stability. Several numerical experiments are presented to support our theoretical results and illustrate the effectiveness of our numerical scheme.