The present study is devoted to investigating the inverse problem of simultaneously reconstructing two source terms that depend solely on time in a system of coupled fractional reaction–diffusion equations. Such coupled systems are fundamental for modelling multispecies anomalous diffusion processes, where the evolution of each state variable is governed by unknown time-varying sources. The inverse problem is tackled using supplementary measurements of the state variables over the spatial domain. The coupling between the two equations presents a distinct complexity, making the simultaneous reconstruction problem notably challenging and important. Results on the unique solvability, are established by employing the Rothe method, whereby the problem is first discretized in time and then stability results for the semi-discrete approximations are derived. These estimates, together with compactness arguments, are then used to rigorously prove the convergence of Rothe approximations to a unique weak solution. Finally, the proposed method's effectiveness and stability are further validated through a series of numerical simulations.