diff --git a/README.md b/README.md
index 614883f..1b5b95a 100644
--- a/README.md
+++ b/README.md
@@ -80,4 +80,72 @@ DEBUG:__main__:
 
 The result is saved in the folder `data/topdata` as a JSON file — giving a model name helps keeping
 track of these outputs.
-Note that for Calabi–Yau dimensions larger than four, the additional 
\ No newline at end of file
+Note that for Calabi–Yau dimensions larger than four, the additional
+
+## period_computation
+
+`period_computation.sage` provides classes for working with Picard–Fuchs operators and their
+period solutions: `PFOperator`, `PFIdeal` and `Period`, together with Ansatz variants of the first and
+last (`PFOperatorAnsatz`, `PeriodAnsatz`) used to search for unknown operators or periods of a given
+z- and theta-degree.
+
+A `PFOperator` is parsed from a string in the variables `z0, ..., z<n-1>` and `theta0, ..., theta<n-1>`,
+the logarithmic derivatives theta_i = z_i d/dz_i. For example, the quintic's Picard&ndash;Fuchs operator:
+
+```python
+L = PFOperator(
+    "theta0^4 - 3125*z0*theta0^4 - 6250*z0*theta0^3 - 4375*z0*theta0^2 - 1250*z0*theta0 - 120*z0",
+    no_variables=1,
+)
+
+L.simplify().operator_string
+```
+
+```term
+'-5*(5*theta0 + 4)*(5*theta0 + 3)*(5*theta0 + 2)*(5*theta0 + 1)*z0 + theta0^4'
+```
+
+An operator (or a `PFIdeal` of several) can be solved for its power series solution at given indicial
+exponents and order.
+
+```python
+ideal = PFIdeal([L])
+period = ideal.find_power_series_solution(indicials=[0], order=3)[0]
+
+period.period_string
+```
+
+```term
+'168168000*z0^3 + 113400*z0^2 + 120*z0 + 1'
+```
+
+The reverse direction is supported too: given a `Period`, `find_annihilating_operators` searches for
+`PFOperator`s of a given z- and theta-degree that annihilate it, by solving an Ansatz of unknown
+coefficients via linear algebra. Both directions extend to several moduli, e.g. for the two-parameter
+model P_{2,2,2,1,1}[8]:
+
+```python
+M1 = PFOperator(
+    "theta1*(-2*theta0 + 2*theta1 - 1) + 2*(theta0 - 2*theta1 - 1)*(theta0 - 2*theta1)*z1",
+    no_variables=2,
+)
+M2 = PFOperator(
+    "theta0^2*(2*(theta0 - 2*theta1)*z1 - theta1) - 16*(2*theta0 + 1)*(4*theta0 + 1)*(4*theta0 + 3)*z0*z1",
+    no_variables=2,
+)
+ideal = PFIdeal([M1, M2])
+
+period = ideal.find_power_series_solution(indicials=[0, 1 / 2], order=6)[0]
+recovered = period.find_annihilating_operators(z_degree=1, theta_degree=2)
+
+period.period_string
+recovered[0].simplify().operator_string
+```
+
+```term
+'-1/45045*(60886425600*z0^3*z1^3 - 2767564800*z0^2*z1^4 + 100638720*z0*z1^5 - 14192640*z1^6 + 830269440*z0^2*z1^3 - 30750720*z0*z1^4 + 4193280*z1^5 - 242161920*z0^2*z1^2 + 9884160*z0*z1^3 - 1281280*z1^4 - 3459456*z0*z1^2 + 411840*z1^3 + 1441440*z0*z1 - 144144*z1^2 + 60060*z1 - 45045)*sqrt(z1)'
+'-2*(theta0 - 2*theta1)*(theta0 - 2*theta1 - 1)*z1 + (2*theta0 - 2*theta1 + 1)*theta1'
+```
+
+`recovered[0]` is, up to scale, `M1` &mdash; recovered purely from `M1`, `M2`'s shared power series
+solution. 
\ No newline at end of file
diff --git a/sage/period_computation.sage b/sage/period_computation.sage
index 56eb473..70d5a65 100644
--- a/sage/period_computation.sage
+++ b/sage/period_computation.sage
@@ -10,7 +10,7 @@ logging.basicConfig(level=logging.DEBUG)
 
 class PFOperator:
     """
-    A class representing a Picard-Fuchs operator in a single variable z.
+    A class representing a Picard-Fuchs operator in a given number of variables (moduli).
     Independent of the given order, the variables are assumed to be left of
     the derivatives.
 
diff --git a/tests/test_period_computation.py b/tests/test_period_computation.py
index 2d368e6..81a2270 100644
--- a/tests/test_period_computation.py
+++ b/tests/test_period_computation.py
@@ -124,7 +124,7 @@ def test_find_power_series_solution_handles_rational_indicial():
     assert ideal.find_power_series_solution(indicials=[0], order=3) == []
 
 
-def test_ideal_finds_holomorphic_solution_and_recovers_first_operator():
+def test_ideal_finds_power_series_solution_and_recovers_first_operator():
     # Picard--Fuchs ideal for P_{22211}[8]:
     #   M1 = theta2*(-2*theta1+2*theta2-1) + 2*(theta1-2*theta2-1)*(theta1-2*theta2)*z2
     #   M2 = theta1^2*(2*(theta1-2*theta2)*z2-theta2) - 16*(2*theta1+1)*(4*theta1+1)*(4*theta1+3)*z1*z2