The Griffiths phase in systems with quenched disorder occurs below the ordering transition of the pure system down to the ordering transition of the actual disordered system. While it does not exhibit long-range order, large fluctuations in the disorder degrees of freedom result in exponentially rare, long-range ordered states and hence the occurrence of broad distributions in response functions. Inside the Griffiths phase of the two-dimensional bond-diluted Ising model the distribution of the magnetic susceptibility is expected to have such a broad, exponential tail. A large-deviation Monte Carlo algorithm is used to sample this distribution and the exponential tail is extracted over a wide range of the support down to very small probabilities of the order of 10^{-300}. We study the behavior of the susceptibility distribution across the full phase diagram, from the paramagnetic state through the Griffiths phase to the ferromagnetically ordered system and down to the zero-temperature point. We extract the rate function of large-deviation theory as well as its finite-size scaling behavior and we reveal interesting differences and similarities between the cases. A connection between the fraction of ferromagnetic bonds in a given disorder sample and the size of the magnetic susceptibility is demonstrated numerically.