diff --git a/CHANGELOG.md b/CHANGELOG.md
index cb702d9df..805e8d22f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,10 +13,11 @@ Changelog
 -   Broadened the scope of `flake8` compliance to the include the `examples` and
     `docs` directories, and thus the whole repository (@tommy-moffat, gh-1113).
 -   `DEFGATE ... AS PAULI-SUM` is now supported (@ecpeterson, gh-1125).
+-   Add unit test for validating Trotterization order (@jmbr, gh-1120).
 
 ### Bugfixes
 
--   Minor fixes for examples/1.3_vqe_demo.py and examples/quantum_walk.ipynb
+-   Minor fixes for `examples/1.3_vqe_demo.py` and `examples/quantum_walk.ipynb`
     (@appleby, gh-1116).
 
 [v2.14](https://github.com/rigetti/pyquil/compare/v2.13.0...v2.14.0) (November 25, 2019)
diff --git a/pyquil/paulis.py b/pyquil/paulis.py
index a49613ea8..e2ec3cf25 100644
--- a/pyquil/paulis.py
+++ b/pyquil/paulis.py
@@ -952,7 +952,7 @@ def suzuki_trotter(trotter_order: int, trotter_steps: int) -> List[Tuple[float,
     to the A operator, o=1 corresponds to the B operator, and w is the
     coefficient in the exponential. For example, a second order Suzuki-Trotter
     approximation to exp(A + B) results in the following
-    [(0.5/trotter_steps, 0), (1/trotteri_steps, 1),
+    [(0.5/trotter_steps, 0), (1/trotter_steps, 1),
     (0.5/trotter_steps, 0)] * trotter_steps.
 
     :param int trotter_order: order of Suzuki-Trotter approximation
diff --git a/pyquil/tests/test_paulis.py b/pyquil/tests/test_paulis.py
index a17ea90c1..e188c522f 100644
--- a/pyquil/tests/test_paulis.py
+++ b/pyquil/tests/test_paulis.py
@@ -29,6 +29,7 @@
                            ID, UnequalLengthWarning, exponentiate, trotterize, is_zero,
                            check_commutation, commuting_sets, term_with_coeff, sI, sX, sY, sZ,
                            ZERO, is_identity)
+from pyquil.unitary_tools import program_unitary
 from pyquil.quil import Program
 
 
@@ -434,6 +435,39 @@ def test_trotterize():
     assert prog == result_prog
 
 
+def test_trotterize_order():
+    def expmi(hermitian_matrix):
+        """Compute the matrix exponential of -1j * hermitian_matrix."""
+        L, Q = np.linalg.eigh(hermitian_matrix)
+        return Q @ np.diag(np.exp(-1j * L)) @ Q.conj().T
+
+    def error(order, time_step_length):
+        a_pauli = time_step_length * sZ(0) * sY(1) * sX(2)
+        a_program = a_pauli.program
+
+        b_pauli = time_step_length * sX(0) * sZ(1) * sY(2)
+        b_program = b_pauli.program
+
+        num_qubits = len(a_program.get_qubits())
+        assert num_qubits == len(b_program.get_qubits())
+
+        a = program_unitary(a_program, num_qubits)
+        b = program_unitary(b_program, num_qubits)
+        a_plus_b = a + b
+        exp_a_plus_b = expmi(time_step_length * a_plus_b)
+
+        trotter_program = trotterize(a_pauli, b_pauli, trotter_order=order)
+        trotter = program_unitary(trotter_program, num_qubits)
+
+        return np.linalg.norm(exp_a_plus_b - trotter, np.inf)
+
+    xs = 10**np.logspace(-1, -6, 10)
+    for order in [1, 2, 3, 4]:
+        ys = [error(order, float(x)) for x in xs]
+        p = np.polyfit(np.log10(xs), np.log10(ys), 1)
+        assert p[0] >= order, f'Bound not satisfied with order={order}: the slope is {p[0]}'
+
+
 def test_is_zero():
     with pytest.raises(TypeError):
         is_zero(1)