We study an inverse problem for a coupled system of semilinear Helmholtz equations where we are interested in reconstructing multiple coefficients in the system from internal data measured in applications such as thermoacoustic imaging. We derive results on the uniqueness and stability of the inverse problem in the case of small boundary data based on the technique of first- and higher-order linearization. Numerical simulations are provided to illustrate the quality of reconstructions that can be expected from noisy data.