  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-deviations 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.
