diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
index 75c9119..f85e5ee 100644
--- a/docs/source/changelog.rst
+++ b/docs/source/changelog.rst
@@ -8,6 +8,7 @@ Unreleased
 
 * `Issue #32 <https://github.com/pmacg/py-sgtl/issues/32>`_ - add tensor product method for combining graphs
 * `Issue #27 <https://github.com/pmacg/py-sgtl/issues/27>`_ - add option to plot spectrum of graph
+* `Issue #31 <https://github.com/pmacg/py-sgtl/issues/31>`_ - add graphs together
 
 0.4.4 - 2022-01-07
 ------------------
diff --git a/docs/source/getting-started.rst b/docs/source/getting-started.rst
index a012ba3..5853c34 100644
--- a/docs/source/getting-started.rst
+++ b/docs/source/getting-started.rst
@@ -50,6 +50,20 @@ Or, you can load a more complex graph from an edgelist file like this.
 See the documentation of the :any:`sgtl.graph.from_edgelist` method for more information on the
 required format of the edgelist file.
 
+You can also add graphs together, which corresponds to adding
+their adjacency matrices.
+
+    >>> import sgtl.graph
+    >>> g1 = sgtl.graph.cycle_graph(5)
+    >>> g2 = sgtl.graph.path_graph(5)
+    >>> g3 = g1 + g2
+    >>> g3.adjacency_matrix().toarray()
+    array([[0., 2., 0., 0., 1.],
+           [2., 0., 2., 0., 0.],
+           [0., 2., 0., 2., 0.],
+           [0., 0., 2., 0., 2.],
+           [1., 0., 0., 2., 0.]])
+
 Viewing the spectrum of a graph
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Investigating the spectrum of a graph is very simple. For example, you could visualise
diff --git a/sgtl/graph.py b/sgtl/graph.py
index a2667f3..1f7c687 100644
--- a/sgtl/graph.py
+++ b/sgtl/graph.py
@@ -66,6 +66,23 @@ def to_networkx(self) -> nx.Graph:
         """
         return nx.Graph(self.adjacency_matrix())
 
+    def __add__(self, other):
+        """
+        Adding two graphs requires that they have the same number of vertices. the sum of the graphs is simply the graph
+        constructed by adding their adjacency matrices together.
+
+        You can also just add a sparse matrix to the graph directly.
+        """
+        if isinstance(other, scipy.sparse.spmatrix):
+            if other.shape[0] != self.number_of_vertices():
+                raise AssertionError("Graphs must have equal number of vertices.")
+            return Graph(self.adjacency_matrix() + other)
+
+        if other.number_of_vertices() != self.number_of_vertices():
+            raise AssertionError("Graphs must have equal number of vertices.")
+
+        return Graph(self.adjacency_matrix() + other.adjacency_matrix())
+
     def draw(self):
         """
         Plot the graph, by first converting to a networkx graph. This will use the default networkx plotting
diff --git a/sgtl/random.py b/sgtl/random.py
index 754046a..dcc94c4 100644
--- a/sgtl/random.py
+++ b/sgtl/random.py
@@ -208,15 +208,16 @@ def ssbm(n: int, k: int, p: float, q: float, directed=False):
     """
     Generate a graph from the symmetric stochastic block model.
 
-    Generates a graph with n vertices and k clusters. Every cluster will have floor(n/k) vertices. The probability of
-    each edge inside a cluster is given by p. The probability of an edge between two different clusters is q.
+    Generates a graph with n vertices and k clusters. Every cluster will have :math:`\\lfloor n/k \\rfloor` vertices.
+    The probability of each edge inside a cluster is given by :math:`p`. The probability of an edge between two
+    different clusters is :math:`q`.
 
     :param n: The number of vertices in the graph.
     :param k: The number of clusters.
     :param p: The probability of an edge inside a cluster.
     :param q: The probability of an edge between clusters.
     :param directed: Whether to generate a directed graph.
-    :return: The generated graph as an ``sgtl.Graph`` object.
+    :return: The generated graph as an :py:class:`sgtl.graph.Graph` object.
     """
     # We are ok with using the 'n', 'k', 'p', and 'q' variable names for the stochastic block model - these are
     # standard notation for this model.
diff --git a/tests/test_graph.py b/tests/test_graph.py
index 6dc0ca0..a05d610 100644
--- a/tests/test_graph.py
+++ b/tests/test_graph.py
@@ -521,3 +521,31 @@ def test_tensor_product():
     adj_mat_diff = new_graph.adjacency_matrix() - expected_adj_mat
     adj_mat_diff.eliminate_zeros()
     assert adj_mat_diff.nnz == 0
+
+
+def test_add():
+    # Test adding two graphs together
+    g1 = sgtl.graph.cycle_graph(5)
+    g2 = sgtl.graph.path_graph(5)
+    g3 = g1 + g2
+
+    expected_adj_mat = sp.sparse.csr_matrix([[0, 2, 0, 0, 1],
+                                             [2, 0, 2, 0, 0],
+                                             [0, 2, 0, 2, 0],
+                                             [0, 0, 2, 0, 2],
+                                             [1, 0, 0, 2, 0]])
+    adj_mat_diff = g3.adjacency_matrix() - expected_adj_mat
+    adj_mat_diff.eliminate_zeros()
+    assert adj_mat_diff.nnz == 0
+
+    # Try adding tha identity matrix
+    g3 = g1 + sp.sparse.eye(5)
+    expected_adj_mat = sp.sparse.csr_matrix([[1, 1, 0, 0, 1],
+                                             [1, 1, 1, 0, 0],
+                                             [0, 1, 1, 1, 0],
+                                             [0, 0, 1, 1, 1],
+                                             [1, 0, 0, 1, 1]])
+    adj_mat_diff = g3.adjacency_matrix() - expected_adj_mat
+    adj_mat_diff.eliminate_zeros()
+    assert adj_mat_diff.nnz == 0
+