diff --git a/docs/gravitational_wave/examples/merger_rates.yml b/docs/gravitational_wave/examples/merger_rates.yml
new file mode 100644
index 000000000..d8d96cc75
--- /dev/null
+++ b/docs/gravitational_wave/examples/merger_rates.yml
@@ -0,0 +1,39 @@
+cosmology: !astropy.cosmology.default_cosmology.get
+z_range: !numpy.arange [0, 2.01, 0.1]
+mag_lim: 30
+sky_area: 1 deg2
+M_star_red: !astropy.modeling.models.Linear1D [-0.70798041, -20.37196157]
+phi_star_red: !astropy.modeling.models.Exponential1D [0.0035097, -1.41649]
+alpha_red: -0.5
+M_star_blue: !astropy.modeling.models.Linear1D [-0.9408582, -20.40492365]
+phi_star_blue: !astropy.modeling.models.Exponential1D [0.00370253, -9.73858]
+alpha_blue: -1.3
+tables:
+  merger_rates:
+    redshift: !skypy.galaxy.redshift.schechter_lf_redshift
+      redshift: $z_range
+      M_star: $M_star_red
+      phi_star: $phi_star_red
+      alpha: $alpha_red
+      m_lim: $mag_lim
+      sky_area: $sky_area
+    magnitude: !skypy.galaxy.luminosity.schechter_lf_magnitude
+      redshift: $merger_rates.redshift
+      M_star: $M_star_red
+      alpha: $alpha_red
+      m_lim: $mag_lim
+    luminosity: !skypy.galaxy.luminosity.luminosity_from_absolute_magnitude
+      absolute_magnitude: $merger_rates.magnitude
+      zeropoint: Lsun_B
+    ns_ns_rate: !skypy.gravitational_wave.b_band_merger_rate
+      population: NS-NS
+      optimism: low
+      luminosity: $merger_rates.luminosity
+    ns_bh_rate: !skypy.gravitational_wave.b_band_merger_rate
+      population: NS-BH
+      optimism: low
+      luminosity: $merger_rates.luminosity
+    bh_bh_rate: !skypy.gravitational_wave.b_band_merger_rate
+      population: BH-BH
+      optimism: low
+      luminosity: $merger_rates.luminosity
diff --git a/docs/gravitational_wave/index.rst b/docs/gravitational_wave/index.rst
index 0bd0f0db1..6e0811a76 100644
--- a/docs/gravitational_wave/index.rst
+++ b/docs/gravitational_wave/index.rst
@@ -2,4 +2,36 @@
 Gravitational Waves (`skypy.gravitational_wave`)
 ************************************************
 
+SkyPy provides methods to model physical processes involving gravitational waves.
+Here we demonstrate calculating the merger rates for three
+different types of compact binary mergers.
+For this we use `~skypy.gravitational_wave.b_band_merger_rate` and
+the configuration file provided in the `~skypy.examples` directory:
+
+
+.. plot::
+   :include-source: false
+
+    from skypy.pipeline import Pipeline
+    pipeline = Pipeline.read('examples/merger_rates.yml')
+    pipeline.execute()
+    merger = pipeline['merger_rates']
+
+    import matplotlib.pyplot as plt
+
+    plt.figure()
+    plt.hist(merger['ns_ns_rate'], histtype='step', bins=np.logspace(-4,2,50), label='NS-NS')
+    plt.hist(merger['ns_bh_rate'], histtype='step', bins=np.logspace(-4,2,50), label='NS-BH')
+    plt.hist(merger['bh_bh_rate'], histtype='step', bins=np.logspace(-4,2,50), label='BH-BH')
+    plt.legend()
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.xlabel('Merger rate, $R\, [year^{-1}$]')
+
+    plt.show()
+
+
+Note: ‘NS-NS’ means neutron star - neutron star, ‘NS-BH’ is neutron star
+- black hole and ‘BH-BH’ is black hole - black hole.
+
 .. automodule:: skypy.gravitational_wave