diff --git a/examples/1INFILTR/D.py b/examples/1INFILTR/D.py
index bbdf49d..037656d 100755
--- a/examples/1INFILTR/D.py
+++ b/examples/1INFILTR/D.py
@@ -6,11 +6,9 @@
 Diffusivity plot.
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 from fronts.D import van_genuchten
-
 from validation import r_unit, t_unit
 
 epsilon = 1e-6
diff --git a/examples/1INFILTR/solve.py b/examples/1INFILTR/solve.py
index 3e756c3..961a5cb 100755
--- a/examples/1INFILTR/solve.py
+++ b/examples/1INFILTR/solve.py
@@ -3,12 +3,10 @@
 """1INFILTR case from Hydrus-1D, horizontal."""
 
 import matplotlib.pyplot as plt
-
+import validation
 from fronts import solve
 from fronts.D import van_genuchten
 
-import validation
-
 epsilon = 1e-6
 
 Ks = 25  # cm/h
diff --git a/examples/1INFILTR/validation.py b/examples/1INFILTR/validation.py
index d05cc78..692b5ed 100755
--- a/examples/1INFILTR/validation.py
+++ b/examples/1INFILTR/validation.py
@@ -7,13 +7,12 @@
 1001 nodes, water content tolerance = 1e-6
 """
 
-import sys
-import os
 import itertools
+import os
+import sys
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 
 r_unit = "cm"
 t_unit = "h"
diff --git a/examples/HF135/D.py b/examples/HF135/D.py
index ddbb90b..1c5d228 100755
--- a/examples/HF135/D.py
+++ b/examples/HF135/D.py
@@ -6,11 +6,9 @@
 Diffusivity plot.
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 from fronts.D import van_genuchten
-
 from validation import r_unit, t_unit
 
 epsilon = 1e-7
diff --git a/examples/HF135/inverse1.py b/examples/HF135/inverse1.py
index 917ce95..3cc423c 100755
--- a/examples/HF135/inverse1.py
+++ b/examples/HF135/inverse1.py
@@ -11,13 +11,11 @@
 Warning: this example takes ~30 seconds to run to completion.
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-
-from fronts import solve, inverse
-from fronts.D import van_genuchten
-
+import numpy as np
 import validation
+from fronts import inverse, solve
+from fronts.D import van_genuchten
 
 epsilon = 1e-7
 
diff --git a/examples/HF135/inverse2.py b/examples/HF135/inverse2.py
index 77ffe4a..06a00a4 100755
--- a/examples/HF135/inverse2.py
+++ b/examples/HF135/inverse2.py
@@ -10,10 +10,8 @@
 """
 
 import matplotlib.pyplot as plt
-
-from fronts import solve, inverse, o
-
 import validation
+from fronts import inverse, o, solve
 
 D_inverse = inverse(o=o(validation.r, validation.t)[::5], samples=validation.theta[::5])
 # Using only a fifth of the points so that it does not run too slow
diff --git a/examples/HF135/radial.py b/examples/HF135/radial.py
index 45ae8f1..07f40ce 100755
--- a/examples/HF135/radial.py
+++ b/examples/HF135/radial.py
@@ -4,14 +4,12 @@
 Radial flow in a HF135 nitrocellulose membrane.
 """
 
-import numpy as np
-import matplotlib.pyplot as plt
-
 from math import pi
 
+import matplotlib.pyplot as plt
+import numpy as np
 from fronts import solve_flowrate
 from fronts.D import van_genuchten
-
 from validation import r_unit, t_unit
 
 epsilon = 1e-7
diff --git a/examples/HF135/refine.py b/examples/HF135/refine.py
index 754c1ef..1553244 100755
--- a/examples/HF135/refine.py
+++ b/examples/HF135/refine.py
@@ -8,12 +8,10 @@
 """
 
 import matplotlib.pyplot as plt
-
+import validation
 from fronts import solve, solve_from_guess
 from fronts.D import van_genuchten
 
-import validation
-
 epsilon = 1e-7
 
 # Wetting of an HF135 membrane, Van Genuchten model
diff --git a/examples/HF135/solve.py b/examples/HF135/solve.py
index 8b56e6e..1897c66 100755
--- a/examples/HF135/solve.py
+++ b/examples/HF135/solve.py
@@ -8,12 +8,10 @@
 """
 
 import matplotlib.pyplot as plt
-
+import validation
 from fronts import solve
 from fronts.D import van_genuchten
 
-import validation
-
 epsilon = 1e-7
 
 # Wetting of an HF135 membrane, Van Genuchten model
diff --git a/examples/HF135/validation.py b/examples/HF135/validation.py
index c6cefd5..9beb3ff 100755
--- a/examples/HF135/validation.py
+++ b/examples/HF135/validation.py
@@ -9,11 +9,11 @@
 
 """
 
-import sys
 import os
+import sys
 
-import numpy as np
 import matplotlib.pyplot as plt
+import numpy as np
 
 _filename = os.path.join(sys.path[0], "groundwaterFoam_results.csv")
 
diff --git a/examples/exact/D.py b/examples/exact/D.py
index cc0703c..32918f1 100755
--- a/examples/exact/D.py
+++ b/examples/exact/D.py
@@ -2,8 +2,8 @@
 
 """Plot of D for the 'exact' validation case."""
 
-import numpy as np
 import matplotlib.pyplot as plt
+import numpy as np
 
 
 def D(theta):
diff --git a/examples/exact/fromguess.py b/examples/exact/fromguess.py
index 9188298..516bb5f 100755
--- a/examples/exact/fromguess.py
+++ b/examples/exact/fromguess.py
@@ -5,9 +5,8 @@
 `fronts.solve_from_guess`) and compares the solutions.
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 from fronts import solve_from_guess
 
 o = np.linspace(0, 20, 100)
diff --git a/examples/exact/solve.py b/examples/exact/solve.py
index 3790064..2f61226 100755
--- a/examples/exact/solve.py
+++ b/examples/exact/solve.py
@@ -5,9 +5,8 @@
 and compares the solutions.
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 from fronts import solve
 
 theta = solve(D="0.5*(1 - log(theta))", i=0, b=1, verbose=2)
diff --git a/examples/powerlaw/D.py b/examples/powerlaw/D.py
index 6fb566d..8584672 100755
--- a/examples/powerlaw/D.py
+++ b/examples/powerlaw/D.py
@@ -4,9 +4,8 @@
 Plot of D in the powerlaw case.
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 from fronts.D import power_law
 
 k = 4.0
diff --git a/examples/powerlaw/inverse.py b/examples/powerlaw/inverse.py
index af1a6b7..b582494 100755
--- a/examples/powerlaw/inverse.py
+++ b/examples/powerlaw/inverse.py
@@ -2,10 +2,9 @@
 
 """Examples of usage of `fronts.solve` and `fronts.inverse`."""
 
-import numpy as np
 import matplotlib.pyplot as plt
-
-from fronts import solve, inverse
+import numpy as np
+from fronts import inverse, solve
 from fronts.D import power_law
 
 k = 4.0
diff --git a/examples/powerlaw/solve.py b/examples/powerlaw/solve.py
index afbc5cf..d086621 100755
--- a/examples/powerlaw/solve.py
+++ b/examples/powerlaw/solve.py
@@ -2,13 +2,11 @@
 
 """Example of usage of `fronts.solve`."""
 
-import numpy as np
 import matplotlib.pyplot as plt
-
+import numpy as np
 from fronts import solve
 from fronts.D import power_law
 
-
 k = 4.0
 ci = 0.1
 cb = 1.0
diff --git a/fronts/D.py b/fronts/D.py
index ff7fe46..329668c 100644
--- a/fronts/D.py
+++ b/fronts/D.py
@@ -1,5 +1,3 @@
-# -*- coding: utf-8 -*-
-
 """D functions."""
 
 import functools
@@ -42,7 +40,6 @@ def constant(D0):
     numerical solvers are necessary. However, it is provided here given that it
     is the simplest supported function.
     """
-
     if D0 <= 0:
         raise ValueError("D0 must be positive")
 
@@ -103,7 +100,6 @@ def from_expr(expr, vectorized=True, max_derivatives=2):
     Users will rarely need to call this function, as all built-in solver
     functions already do so themselves when they receive an expression as `D`.
     """
-
     expr = sympy.sympify(expr)
 
     free = expr.free_symbols
@@ -181,7 +177,7 @@ def D(theta, derivatives=0):
                 return f01(theta)
 
             if derivatives == 2 and max_derivatives == 2:
-                return f01(theta) + [f2(theta)]
+                return (*f01(theta), f2(theta))
 
             raise ValueError(
                 f"derivatives must be one of {{{', '.join(str(n) for n in range(max_derivatives+1))}}}"
@@ -251,7 +247,9 @@ def D(theta, derivatives=0):
 
 def _as_Ks(Ks=None, k=None, nu=1e-6, g=9.81):
     r"""
-    Return the saturated hydraulic conductivity, computed from the instrinsic
+    Return the saturated hydraulic conductivity.
+
+    The saturated hydraulic conductivity is computed from the instrinsic
     permeability if necessary.
 
     Parameters
@@ -372,10 +370,9 @@ def brooks_and_corey(
 
     References
     ----------
-    [1] BROOKS, R.; COREY, T. Hydraulic properties of porous media. Hydrology
+    [1] BROOKS, R.; COREY, T. Hydraulic properties of porous media. Hydrology
     Papers, Colorado State University, 1964, vol. 24, p. 37.
     """
-
     if alpha <= 0:
         raise ValueError("alpha must be positive")
 
@@ -496,11 +493,10 @@ def van_genuchten(
 
     References
     ----------
-    [1] VAN GENUCHTEN, M. Th. A closed-form equation for predicting the
+    [1] VAN GENUCHTEN, M. Th. A closed-form equation for predicting the
     hydraulic conductivity of unsaturated soils. Soil Science Society of
     America Journal, 1980, vol. 44, no 5, p. 892-898.
     """
-
     if n is not None:
         if m is not None:
             raise TypeError("cannot pass both n and m")
@@ -605,6 +601,8 @@ def letxs(
     theta_range=(0.0, 1.0),
 ):
     r"""
+    Return a LETx + LETs diffusivity function.
+
     Return a diffusivity function that combines the LETx relative permeability
     correlation and the LETs capillary pressure correlation for spontaneous
     imbibition. Both correlations are part of the LET family of hydraulic
@@ -699,7 +697,6 @@ def letxs(
     capillary imbibition models in paper-based microfluidic applications.
     Transport in Porous Media, 2022, vol. 141, no. 7, pp. 1-20.
     """
-
     Ks = _as_Ks(Ks=Ks, k=k, nu=nu, g=g)
 
     # - Code generated with functionstr() from ../symbolic/generate.py - #
@@ -997,7 +994,6 @@ def richards(C, kr, Ks=None, k=None, nu=1e-6, g=9.81):
         In all cases, the argument ``theta`` may be a single float or a NumPy
         array.
     """
-
     Ks = _as_Ks(Ks=Ks, k=k, nu=nu, g=g)
 
     def D(theta, derivatives=0):
diff --git a/fronts/__init__.py b/fronts/__init__.py
index c254fd0..bda149d 100644
--- a/fronts/__init__.py
+++ b/fronts/__init__.py
@@ -1,13 +1,10 @@
-"""
-Numerical library for nonlinear diffusion problems based on the Boltzmann
-transformation.
-"""
+"""Numerical library for nonlinear diffusion problems based on the Boltzmann transformation."""
 
 __version__ = "1.2.2"
 
-from ._boltzmann import ode, BaseSolution, o, do_dr, do_dt, r, t, as_o
-from ._semiinfinite import solve, solve_flowrate, solve_from_guess, Solution
+from ._boltzmann import BaseSolution, as_o, do_dr, do_dt, o, ode, r, t
 from ._inverse import inverse, sorptivity
+from ._semiinfinite import Solution, solve, solve_flowrate, solve_from_guess
 
 __all__ = [
     "solve",
diff --git a/fronts/_boltzmann.py b/fronts/_boltzmann.py
index 6e5e59d..0c5af95 100644
--- a/fronts/_boltzmann.py
+++ b/fronts/_boltzmann.py
@@ -1,4 +1,6 @@
 """
+Internal Boltzmann transformation module.
+
 This module defines the Boltzmann variable, transformation of the partial
 differential equation into the ODE in terms of that variable, and inverse
 transformation of the ODE's solution into a solution to the partial
@@ -32,7 +34,7 @@ def o(r, t):
         The return is a float if both `r` and `t` are floats. Otherwise it is
         an array of the shape that results from broadcasting `r` and `t`.
 
-    See also
+    See Also
     --------
     do_dr
     do_dt
@@ -64,7 +66,7 @@ def do_dr(r, t):
         The return is a float if both `r` and `t` are floats. Otherwise it is
         an array of the shape that results from broadcasting `r` and `t`.
 
-    See also
+    See Also
     --------
     o
     do_dt
@@ -93,7 +95,7 @@ def do_dt(r, t):
         The return is a float if both `r` and `t` are floats. Otherwise it is
         an array of the shape that results from broadcasting `r` and `t`.
 
-    See also
+    See Also
     --------
     o
     do_dr
@@ -120,7 +122,7 @@ def r(o, t):
         The return is a float if both `o` and `t` are floats. Otherwise it is
         an array of the shape that results from broadcasting `o` and `t`.
 
-    See also
+    See Also
     --------
     o
     t
@@ -147,7 +149,7 @@ def t(o, r):
         The return is a float if both `o` and `r` are floats. Otherwise it is
         an array of the shape that results from broadcasting `o` and `r`.
 
-    See also
+    See Also
     --------
     o
     r
@@ -160,9 +162,10 @@ def t(o, r):
 
 def as_o(r=None, t=None, o=None):
     """
-    Transform to the Boltzmann variable if called with `r` and `t`. Passes the
-    values through if called with `o` only. On other combinations of arguments,
-    it raises a `TypeError` with a message explaining valid usage.
+    Transform to the Boltzmann variable if called with `r` and `t`.
+
+    Passes the values through if called with `o` only. On other combinations of
+    arguments, it raises a `TypeError` with a message explaining valid usage.
 
     This function is a helper to define other functions that may be called
     either with `r` and `t`, or with just `o`.
@@ -185,7 +188,7 @@ def as_o(r=None, t=None, o=None):
         Passes `o` through if it is given. Otherwise, it calls the function
         :func:`o` and returns ``o(r,t)``.
 
-    See also
+    See Also
     --------
     o
     """
@@ -267,7 +270,7 @@ def ode(D, radial=False, catch_errors=False):
         and (2,n) respectively, and the return is a NumPy array of shape
         (2,2,n).
 
-    Other parameters
+    Other Parameters
     ----------------
     catch_errors : bool, optional
         Whether to catch exceptions that may be attributed to a domain error of
@@ -296,7 +299,7 @@ def ode(D, radial=False, catch_errors=False):
         If False (default), the exceptions will be allowed to propagate to the
         callers of `fun` and `jac`.
 
-    See also
+    See Also
     --------
     BaseSolution
     o
@@ -306,7 +309,6 @@ def ode(D, radial=False, catch_errors=False):
     If `radial` is other than `False`, the PDE is undefined at :math:`r=0`, and
     therefore the returned ODE is also undefined for :math:`o=0`.
     """
-
     try:
         k = _k[radial]
     except KeyError:
@@ -352,7 +354,7 @@ def fun(o, y):
     def jac(o, y):
         theta, dtheta_do = y
 
-        J = np.empty((2, 2) + np.shape(o))
+        J = np.empty((2, 2, *np.shape(o)))
         J[0, 0] = 0
         J[0, 1] = 1
 
@@ -430,7 +432,7 @@ class BaseSolution:
         Function to evaluate :math:`D` at arbitrary values of the solution.
         Must be callable with a float or NumPy array as its argument.
 
-    See also
+    See Also
     --------
     ode
     """
diff --git a/fronts/_inverse.py b/fronts/_inverse.py
index ab12e39..9d33cc3 100644
--- a/fronts/_inverse.py
+++ b/fronts/_inverse.py
@@ -50,7 +50,7 @@ def inverse(o, samples):
         :math:`D` is guaranteed to be continuous; however, its derivatives are
         not.
 
-    See also
+    See Also
     --------
     o
 
@@ -70,7 +70,6 @@ def inverse(o, samples):
     [2] BRUCE, R. R.; KLUTE, A. The measurement of soil moisture diffusivity.
     Soil Science Society of America Journal, 1956, vol. 20, no. 4, pp. 458-462.
     """
-
     o = np.asarray(o)
 
     if not np.all(np.diff(o) > 0):
diff --git a/fronts/_rootfinding.py b/fronts/_rootfinding.py
index 164317e..6505057 100644
--- a/fronts/_rootfinding.py
+++ b/fronts/_rootfinding.py
@@ -41,6 +41,8 @@ def __init__(
 
 class IterationLimitReached(RuntimeError):
     """
+    Iteration limit reached during a call to a function in this module.
+
     Exception raised when a function in this module does not finish within the
     specified maximum number of iterations.
 
@@ -66,8 +68,7 @@ def bracket_root(
     f, interval, growth_factor=2, maxiter=100, f_interval=(None, None), ftol=None
 ):
     """
-    Find an interval that brackets a root of a function by searching in one
-    direction.
+    Find an interval that brackets a root of a function.
 
     Starting from an interval, it moves and expands the interval in the
     direction of the second endpoint until the interval brackets a root of the
@@ -106,7 +107,7 @@ def bracket_root(
         the result will also include a bracket only if one was found at the
         same time as the root.
 
-    See also
+    See Also
     --------
     bisect
 
@@ -198,6 +199,8 @@ def bracket_root(
 
 class NotABracketError(ValueError):
     """
+    Exception raised when a bracket is not valid.
+
     Exception raised by :func:`bisect` when the interval passed as `bracket`
     does not actually contain a root.
 
@@ -253,7 +256,7 @@ def bisect(f, bracket, ftol=1e-12, maxiter=100, f_bracket=(None, None)):
     result : Result
         Contains the root and the final bracket.
 
-    See also
+    See Also
     --------
     bracket_root : Search for a bracket.
 
@@ -265,7 +268,6 @@ def bisect(f, bracket, ftol=1e-12, maxiter=100, f_bracket=(None, None)):
     qualify as roots, the one where the absolute value of `f` is lower is
     returned.
     """
-
     if ftol < 0:
         raise ValueError("ftol cannot be negative")
 
diff --git a/fronts/_semiinfinite.py b/fronts/_semiinfinite.py
index 8a541d8..06fbdeb 100644
--- a/fronts/_semiinfinite.py
+++ b/fronts/_semiinfinite.py
@@ -1,4 +1,6 @@
 """
+Semi-infinite domain problems with the Boltzmann transformation.
+
 This module uses the Boltzmann transformation to deal with initial-boundary
 value problems in semi-infinite domains.
 """
@@ -7,11 +9,11 @@
 from time import process_time
 
 import numpy as np
-from scipy.integrate import solve_ivp, solve_bvp
+from scipy.integrate import solve_bvp, solve_ivp
 
-from ._boltzmann import ode, BaseSolution, r
-from ._rootfinding import bracket_root, bisect, NotABracketError
-from .D import from_expr, _checked
+from ._boltzmann import BaseSolution, ode, r
+from ._rootfinding import NotABracketError, bisect, bracket_root
+from .D import _checked, from_expr
 
 
 class Solution(BaseSolution):
@@ -79,8 +81,12 @@ def ob(self):
 
     @property
     def oi(self):
-        """float: Value of the Boltzmann variable at which the solution can be
-        considered to be equal to the initial condition."""
+        """
+        float: Parameter :math:`o_i`.
+
+        The value of the Boltzmann variable at which the solution can be
+        considered to be equal to the initial condition.
+        """
         return self._oi
 
     def rb(self, t):
@@ -167,7 +173,9 @@ def fluxb(self, t):
     @property
     def d_dob(self):
         """
-        float: Derivative of the solution with respect to the Boltzmann
+        float: boundary Boltzmann derivative.
+
+        Derivative of the solution with respect to the Boltzmann
         variable at the boundary.
         """
         return self.d_do(o=self.ob)
@@ -227,9 +235,7 @@ class _Shooter:
 
     @staticmethod
     def _native_float_inputs(f):
-        """
-        Speeds up arithmetic by converting NumPy inputs to native floats.
-        """
+        """Speeds up arithmetic by converting NumPy inputs to native floats."""
 
         def wrapper(o, y):
             return f(float(o), (float(y[0]), float(y[1])))
@@ -351,14 +357,17 @@ def integrate(self, b, d_dob):
 
     class ShotLimitReached(RuntimeError):
         """
+        Shot limit reached.
+
         Exception raised when `shoot` is called after the maximum number of
         shots has been reached.
         """
 
     def shoot(self, *args, **kwargs):
         """
-        Calls `integrate` and returns the result's `i_residual`. Each call
-        increments the number of shots.
+        Call `integrate` and returns the result's `i_residual`.
+
+        Each call increments the number of shots.
 
         It raises a `ShotLimitReached` exception if the maximum number of
         allowed shots has been reached.
@@ -536,7 +545,7 @@ def solve(
                 the solver; in particular, this field is never `None` if a
                 `d_dob_bracket` is passed when calling the function.
 
-    Other parameters
+    Other Parameters
     ----------------
     itol : float, optional
         Absolute tolerance for the initial condition.
@@ -574,7 +583,7 @@ def solve(
             * 1: display a termination report
             * 2: also display progress during iterations
 
-    See also
+    See Also
     --------
     solve_from_guess
     solve_flowrate
@@ -910,7 +919,7 @@ def solve_flowrate(
                 solver; in particular, this field is never `None` if a
                 `b_bracket` is passed when calling the function.
 
-    Other parameters
+    Other Parameters
     ----------------
     itol : float, optional
         Absolute tolerance for the initial condition.
@@ -946,7 +955,7 @@ def solve_flowrate(
             * 1: display a termination report
             * 2: also display progress during iterations
 
-    See also
+    See Also
     --------
     solve
 
@@ -1229,7 +1238,7 @@ def solve_from_guess(D, i, b, o_guess, guess, radial=False, max_nodes=1000, verb
             *   `rms_residuals` *(numpy.ndarray, shape (n-1,))* --
                 RMS values of the relative residuals over each mesh interval.
 
-    Other parameters
+    Other Parameters
     ----------------
     max_nodes : int, optional
         Maximum allowed number of mesh nodes.
@@ -1240,7 +1249,7 @@ def solve_from_guess(D, i, b, o_guess, guess, radial=False, max_nodes=1000, verb
             * 1: display a termination report
             * 2: also display progress during iterations
 
-    See also
+    See Also
     --------
     solve
 
diff --git a/symbolic/generate.py b/symbolic/generate.py
index 6fb58a1..ee1c71a 100644
--- a/symbolic/generate.py
+++ b/symbolic/generate.py
@@ -1,7 +1,7 @@
-import sympy
-
 import itertools
 
+import sympy
+
 
 def _derivative_names(var):
     yield "D"
diff --git a/symbolic/letd.py b/symbolic/letd.py
index d63774d..7387490 100755
--- a/symbolic/letd.py
+++ b/symbolic/letd.py
@@ -24,6 +24,6 @@
 Swp = (theta - theta_range[0]) / (theta_range[1] - theta_range[0])
 D = Dwt * Swp**L / (Swp**L + E * (1 - Swp) ** T)
 
-print("D={}".format(D))
+print(f"D={D}")
 
 print(functionstr(theta, D.simplify()))
diff --git a/symbolic/letxs.py b/symbolic/letxs.py
index 10e8d70..8e155a9 100755
--- a/symbolic/letxs.py
+++ b/symbolic/letxs.py
@@ -48,6 +48,6 @@
 
 ################################
 
-print("D={}".format(D))
+print(f"D={D}")
 
 print(functionstr(theta, D))
diff --git a/tests/test_D.py b/tests/test_D.py
index fe7e19b..3d1336d 100644
--- a/tests/test_D.py
+++ b/tests/test_D.py
@@ -1,12 +1,10 @@
-import pytest
-
 import functools
 
+import fronts.D
 import numpy as np
+import pytest
 from autograd import deriv
 
-import fronts.D
-
 
 def van_genuchten_D(theta, m, l, alpha, Ks, theta_range):  # noqa: E741
     # Reference: Van Genuchten (1980) Equation 11
diff --git a/tests/test_examples.py b/tests/test_examples.py
index b9a3220..f39462f 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -1,7 +1,7 @@
-import pytest
-
 import os
 
+import pytest
+
 examples = []
 for root, dirs, files in os.walk("examples"):
     for file in files:
diff --git a/tests/test_inverse.py b/tests/test_inverse.py
index 7361e9e..1db581f 100644
--- a/tests/test_inverse.py
+++ b/tests/test_inverse.py
@@ -1,10 +1,8 @@
-import pytest
-
+import fronts
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
-import fronts
-
 
 def test_exact():
     # Reference: Philip (1960) Table 1, No. 13
diff --git a/tests/test_rootfinding/test_bisect.py b/tests/test_rootfinding/test_bisect.py
index 0661162..328729b 100644
--- a/tests/test_rootfinding/test_bisect.py
+++ b/tests/test_rootfinding/test_bisect.py
@@ -1,10 +1,8 @@
-import pytest
-
 import sys
 
 import fronts._rootfinding as rootfinding
-
-from checkobj import check_result, check_notabracketerror, check_iterationlimitreached
+import pytest
+from checkobj import check_iterationlimitreached, check_notabracketerror, check_result
 
 
 def f(x):
diff --git a/tests/test_rootfinding/test_bracket_root.py b/tests/test_rootfinding/test_bracket_root.py
index ac9ee77..054a15b 100644
--- a/tests/test_rootfinding/test_bracket_root.py
+++ b/tests/test_rootfinding/test_bracket_root.py
@@ -1,9 +1,6 @@
-import pytest
-
-
 import fronts._rootfinding as rootfinding
-
-from checkobj import check_result, check_iterationlimitreached
+import pytest
+from checkobj import check_iterationlimitreached, check_result
 
 
 def f(x):
diff --git a/tests/test_solve.py b/tests/test_solve.py
index 4aaacc5..d7b5ce8 100644
--- a/tests/test_solve.py
+++ b/tests/test_solve.py
@@ -1,10 +1,8 @@
-import pytest
-
-import numpy as np
-from numpy.testing import assert_allclose
-
 import fronts
 import fronts.D
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
 
 
 def test_nogradient():
diff --git a/tests/test_solve_flowrate.py b/tests/test_solve_flowrate.py
index 3df16ef..4d75d43 100644
--- a/tests/test_solve_flowrate.py
+++ b/tests/test_solve_flowrate.py
@@ -1,10 +1,8 @@
-import pytest
-
+import fronts
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
-import fronts
-
 
 def test_Qb():
     Qb = 2
diff --git a/tests/test_solve_from_guess.py b/tests/test_solve_from_guess.py
index 04c7181..e3a0885 100644
--- a/tests/test_solve_from_guess.py
+++ b/tests/test_solve_from_guess.py
@@ -1,8 +1,7 @@
+import fronts
 import numpy as np
 from numpy.testing import assert_allclose
 
-import fronts
-
 
 def test_exact():
     # Reference: Philip (1960) Table 1, No. 13