diff --git a/doc/api.rst b/doc/api.rst
index 19f4e14d8ee242a956fb3a69134c00fd7999c303..85c89e3e9f2461a63be14e3b0865221b8aedb517 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -16,11 +16,16 @@ API
    :members:
    :undoc-members:
 
+.. autofunction:: usadelndsoc.bcs.BCS_Delta
+
 .. automodule:: usadelndsoc
    :members:
    :undoc-members:
 
 .. py:data:: usadelndsoc.logger
 
-   `logging.Logger` instance used.
+   `logging` instance used.
 
+.. automodule:: usadelndsoc.plotting
+   :members:
+   :undoc-members:
diff --git a/usadelndsoc/__init__.py b/usadelndsoc/__init__.py
index de55f55444045f8c2ee9b91c736a703579e4a7b2..1459a04734ab1548b6cbe40687c9cd35673e872b 100644
--- a/usadelndsoc/__init__.py
+++ b/usadelndsoc/__init__.py
@@ -37,6 +37,11 @@ def _init_logging():
 
 
 class with_log_level(contextlib.ContextDecorator):
+    """
+    Context decorator that temporarily changes
+    the log level of `usadelndsoc.logger`.
+    """
+
     def __init__(self, level):
         self._level = level
 
diff --git a/usadelndsoc/bcs.py b/usadelndsoc/bcs.py
index c3d4ffb68e8a68b54ee615dad9d177e0d0f4a464..2410563010d2dad4251b754aef3868907cdc8fb7 100644
--- a/usadelndsoc/bcs.py
+++ b/usadelndsoc/bcs.py
@@ -18,6 +18,20 @@ __all__ = ["BCS_Delta"]
 def BCS_Delta(T, Tc, h=0):
     """
     Evaluate Delta in a bulk BCS superconductor.
+
+    Parameters
+    ----------
+    T : float
+        Temperature
+    Tc : float
+        BCS critical temperature
+    h : float
+        Exchange field
+
+    Returns
+    -------
+    Delta : float
+        Order parameter
     """
     Delta_0 = pi / exp(euler_gamma) * Tc
 
diff --git a/usadelndsoc/plotting.py b/usadelndsoc/plotting.py
index 93bf6b70cc18fb50e99545e24149ffdd54aeac26..e8e3e720d189d9468b2e8346fa6312a50b4bc088 100644
--- a/usadelndsoc/plotting.py
+++ b/usadelndsoc/plotting.py
@@ -17,6 +17,9 @@ __all__ = [
 
 
 def plot_Delta_J_2d(x, y, Delta, J, ax=None):
+    """
+    Color plot of order parameter and supercurrents.
+    """
     if ax is None:
         ax = plt.gca()
 
@@ -40,6 +43,10 @@ def plot_Delta_J_2d(x, y, Delta, J, ax=None):
 
 
 class ComplexColorMesh:
+    """
+    Object representing Matplotlib complex colormap.
+    """
+
     def __init__(self, phase, amplitude):
         self.phase = phase
         self.amplitude = amplitude
@@ -62,6 +69,22 @@ class ComplexColorMesh:
 
 
 def pcolormesh_complex(x, y, data, ax=None):
+    """
+    Plot complex-valued function as a color map.
+
+    Lightness indicates amplitude, color indicates phase angle.
+
+    Parameters
+    ----------
+    x : array of float, shape (nx,)
+    y : array of float, shape (ny,)
+    data : array of complex, shape (nx, ny)
+    ax : matplotlib axis
+
+    Returns
+    -------
+    cmesh : ComplexColorMesh
+    """
     nx = data.shape[0]
     ny = data.shape[1]
     dx = np.median(np.diff(x))
@@ -98,6 +121,7 @@ def pcolormesh_complex(x, y, data, ax=None):
 
 @functools.lru_cache()
 def cmap_hue_cyclic():
+    """Cyclic hue Matplotlib color map"""
     from ._plotting_data import data
 
     return _get_cmap(data)
@@ -105,6 +129,7 @@ def cmap_hue_cyclic():
 
 @functools.lru_cache()
 def cmap_huegray_r():
+    """Reverse grayscale colormap, corresponding to cmap_hue_cyclic()"""
     from ._plotting_data import data_l
 
     return _get_cmap(np.c_[data_l, data_l, data_l][::-1])
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index 2da894a72d3464c5b5e2a3c7ed0388a082d332bd..86883303a68bb02bb0082267d86054cb2442b034 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -1200,6 +1200,31 @@ def cpr(
     filename=None,
     **selfcons_kw,
 ):
+    """
+    Compute current-phase relation with selfconsistency.
+
+    Adjust superconducting phases in cells indicated by `phase_mask`.
+
+    Parameters
+    ----------
+    solve : Solver
+        Solver instance to use
+    T : float
+        Temperature
+    T_c0 : array of float
+        Critical temperature
+    phase_mask : array of bool
+        Cells where to tune superconducting phase
+    max_points : int
+    ds : float
+        Step size
+    phi0 : float
+    filename : str, optional
+        File to resume or save progress. Computation is resumed if
+        the file exists and parameter values in it match.
+    **selfcons_kw
+        Parameters passed on to `Solver.self_consistency`.
+    """
     nx, ny = solver.shape
 
     phase_mask = phase_mask.astype(bool)
diff --git a/usadelndsoc/util.py b/usadelndsoc/util.py
index 15766c13de62029d893def37e63596888f855779..f256b08754e0a7c4ce28acc4e3de7d973a70c782 100644
--- a/usadelndsoc/util.py
+++ b/usadelndsoc/util.py
@@ -83,6 +83,28 @@ def vectorize_parallel(
     excluded=None,
     noarray=False,
 ):
+    """
+    Decorator to vectorize and function to evaluate elements in parallel.
+    Uses *joblib* for parallelization.
+
+    Parameters
+    ----------
+    func : function
+    backend : str
+        Joblib backend
+    n_jobs : int
+        Number of jobs, -1 means number of CPUs.
+    returns_object : bool
+        Put function return values to an object array.
+    batch_size : {"auto", int}
+    included : tuple
+        Function argument names to include in vectorization.
+    excluded : tuple
+        Function argument names to exclude from vectorization.
+    noarray : bool
+        Convert array-valued return values to scalars.
+    """
+
     if func is None:
 
         def deco(func):