diff --git a/doc/conf.py b/doc/conf.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_ZG9jL2NvbmYucHk=..c7d429d47f2e44566ad76e30b8977572b074db44_ZG9jL2NvbmYucHk= 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -21,7 +21,8 @@
 
 from fluiddoc import mock_modules
 
-mock_modules(('scipy', 'scipy.sparse', 'scipy.sparse.linalg', 'scipy.fftpack'))
+mock_modules(('scipy', 'scipy.sparse', 'scipy.sparse.linalg', 'scipy.fftpack',
+              'basilisk', 'basilisk.stream',))
 
 from fluiddoc.ipynb_maker import ipynb_to_rst
 ipynb_to_rst()
diff --git a/doc/examples.rst b/doc/examples.rst
index aa00499a506c22ee9d6698e90b106bb5afec06bd_ZG9jL2V4YW1wbGVzLnJzdA==..c7d429d47f2e44566ad76e30b8977572b074db44_ZG9jL2V4YW1wbGVzLnJzdA== 100644
--- a/doc/examples.rst
+++ b/doc/examples.rst
@@ -11,3 +11,4 @@
    examples/init-fields-in-script-3d
    examples/forcing-in-script
    examples/static-layers
+   examples/running-third-party-basilisk
diff --git a/doc/examples/running-third-party-basilisk.rst b/doc/examples/running-third-party-basilisk.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c7d429d47f2e44566ad76e30b8977572b074db44_ZG9jL2V4YW1wbGVzL3J1bm5pbmctdGhpcmQtcGFydHktYmFzaWxpc2sucnN0
--- /dev/null
+++ b/doc/examples/running-third-party-basilisk.rst
@@ -0,0 +1,10 @@
+Run a third party solver using FluidSim: Basilisk
+=================================================
+
+A proof of concept on how FluidSim can act as a framework to extend third party
+solvers. Provided `Basilisk <basilisk.fr>`__ is installed and available for import
+the script can be launched as ``python simul_third_party_basilisk.py`` from the
+``doc/examples`` directory.
+
+.. literalinclude:: simul_third_party_basilisk.py
+
diff --git a/fluidsim/base/basilisk/run_example.py b/doc/examples/simul_third_party_basilisk.py
similarity index 100%
rename from fluidsim/base/basilisk/run_example.py
rename to doc/examples/simul_third_party_basilisk.py
diff --git a/fluidsim/__init__.py b/fluidsim/__init__.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vX19pbml0X18ucHk=..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vX19pbml0X18ucHk= 100644
--- a/fluidsim/__init__.py
+++ b/fluidsim/__init__.py
@@ -20,6 +20,7 @@
    operators
    solvers
    magic
+   pkgload
 
 """
 
@@ -33,6 +34,7 @@
 path_dir_results = Path(FLUIDSIM_PATH)
 
 from .util.util import (
+    available_solver_keys,
     import_module_solver_from_key,
     import_simul_class_from_key,
     load_sim_for_plot,
@@ -44,6 +46,7 @@
 
 from .base.params import load_params_simul
 
+
 # clean up
 from . import util
 
@@ -47,6 +50,7 @@
 # clean up
 from . import util
 
+
 del util
 
 
diff --git a/fluidsim/base/__init__.py b/fluidsim/base/__init__.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9fX2luaXRfXy5weQ==..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/base/__init__.py
+++ b/fluidsim/base/__init__.py
@@ -16,5 +16,13 @@
    output
    forcing
    preprocess
+   sphericalharmo
+
+Third party solver interfaces:
+
+.. autosummary::
+   :toctree:
+
+   basilisk
 
 """
diff --git a/fluidsim/base/basilisk/__init__.py b/fluidsim/base/basilisk/__init__.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9fX2luaXRfXy5weQ==..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/base/basilisk/__init__.py
+++ b/fluidsim/base/basilisk/__init__.py
@@ -0,0 +1,17 @@
+"""Basilisk solvers (:mod:`fluidsim.base.basilisk`)
+===================================================
+A reference implementation to use 2D solvers in `Basilisk <basilisk.fr>`__ via
+FluidSim.
+
+Provides:
+
+.. autosummary::
+   :toctree:
+
+    solver
+    state
+    operators
+    time_stepping
+    output
+
+"""
diff --git a/fluidsim/base/basilisk/operators.py b/fluidsim/base/basilisk/operators.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vcGVyYXRvcnMucHk=..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vcGVyYXRvcnMucHk= 100644
--- a/fluidsim/base/basilisk/operators.py
+++ b/fluidsim/base/basilisk/operators.py
@@ -1,3 +1,13 @@
+"""Operators for Basilisk (:mod:`fluidsim.base.basilisk.operators`)
+===================================================================
+
+Provides:
+
+.. autoclass:: OperatorsBasilisk2D
+   :members:
+   :private-members:
+
+"""
 from __future__ import division
 
 from builtins import object
diff --git a/fluidsim/base/basilisk/output/__init__.py b/fluidsim/base/basilisk/output/__init__.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vdXRwdXQvX19pbml0X18ucHk=..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9vdXRwdXQvX19pbml0X18ucHk= 100644
--- a/fluidsim/base/basilisk/output/__init__.py
+++ b/fluidsim/base/basilisk/output/__init__.py
@@ -1,3 +1,14 @@
+"""Module to generate output from Basilisk (:mod:`fluidsim.base.basilisk.output`)
+==================================================================================
+
+Provides:
+
+.. autoclass:: OutputBasilisk
+   :members:
+   :private-members:
+
+"""
+
 from ...output.base import OutputBase  # , SpecificOutput
 
 
diff --git a/fluidsim/base/basilisk/solver.py b/fluidsim/base/basilisk/solver.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9zb2x2ZXIucHk=..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9zb2x2ZXIucHk= 100644
--- a/fluidsim/base/basilisk/solver.py
+++ b/fluidsim/base/basilisk/solver.py
@@ -1,3 +1,17 @@
+"""Basilisk solver (:mod:`fluidsim.base.basilisk.solver`)
+=========================================================
+
+Provides:
+
+.. autoclass:: InfoSolverBasilisk
+   :members:
+   :private-members:
+
+.. autoclass:: SimulBasilisk
+   :members:
+   :private-members:
+
+"""
 from __future__ import print_function
 
 from fluidsim.base.solvers.base import SimulBase
@@ -7,6 +21,7 @@
 
 
 class InfoSolverBasilisk(InfoSolverBase):
+    """Contain the information on a Basilisk solver."""
     def _init_root(self):
 
         super(InfoSolverBasilisk, self)._init_root()
@@ -33,6 +48,7 @@
 
 
 class SimulBasilisk(SimulBase):
+    """A solver for Basilisk."""
     InfoSolver = InfoSolverBasilisk
 
     def __init__(self, params):
@@ -36,6 +52,7 @@
     InfoSolver = InfoSolverBasilisk
 
     def __init__(self, params):
+        """Initialize parameters, state fields, and event loop."""
         bas = self.basilisk = basilisk
         super(SimulBasilisk, self).__init__(params)
 
diff --git a/fluidsim/base/basilisk/state.py b/fluidsim/base/basilisk/state.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9zdGF0ZS5weQ==..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay9zdGF0ZS5weQ== 100644
--- a/fluidsim/base/basilisk/state.py
+++ b/fluidsim/base/basilisk/state.py
@@ -1,3 +1,15 @@
+"""State of the variables (:mod:`fluidsim.base.basilisk.state`)
+===============================================================
+
+Provides:
+
+.. autoclass:: StateBasilisk
+   :members:
+   :private-members:
+
+"""
+
+
 from ..state import StateBase
 
 
@@ -17,7 +29,7 @@
         )
 
     def compute(self, key):
-
+        """Compute scalar fields such a component of the velocity or vorticity."""
         if key == "ux":
             scalar = self.sim.basilisk.uf.x
         elif key == "uy":
diff --git a/fluidsim/base/basilisk/time_stepping.py b/fluidsim/base/basilisk/time_stepping.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay90aW1lX3N0ZXBwaW5nLnB5..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9iYXNpbGlzay90aW1lX3N0ZXBwaW5nLnB5 100644
--- a/fluidsim/base/basilisk/time_stepping.py
+++ b/fluidsim/base/basilisk/time_stepping.py
@@ -1,3 +1,14 @@
+"""Time stepping (:mod:`fluidsim.base.basilisk.time_stepping`)
+==============================================================
+
+
+Provides:
+
+.. autoclass:: TimeSteppingBasilisk
+   :members:
+   :private-members:
+
+"""
 from __future__ import print_function
 
 from time import time
@@ -7,6 +18,9 @@
 
 
 class TimeSteppingBasilisk(object):
+    """Time stepping class to handle Basilisk's event loop and FluidSim output.
+
+    """
     @staticmethod
     def _complete_params_with_default(params):
         """This static method is used to complete the *params* container.
diff --git a/fluidsim/base/sphericalharmo/solver.py b/fluidsim/base/sphericalharmo/solver.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9zb2x2ZXIucHk=..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vYmFzZS9zcGhlcmljYWxoYXJtby9zb2x2ZXIucHk= 100644
--- a/fluidsim/base/sphericalharmo/solver.py
+++ b/fluidsim/base/sphericalharmo/solver.py
@@ -1,3 +1,3 @@
-"""Base solver (:mod:`fluidsim.base.solvers.spherical_harmo`)
-=============================================================
+"""Base solver (:mod:`fluidsim.base.solvers.sphericalharmo`)
+============================================================
 
@@ -3,5 +3,5 @@
 
-.. autoclass:: InfoSolver
+.. autoclass:: InfoSolverSphericalHarmo
    :members:
    :private-members:
 
diff --git a/fluidsim/solvers/__init__.py b/fluidsim/solvers/__init__.py
index aa00499a506c22ee9d6698e90b106bb5afec06bd_Zmx1aWRzaW0vc29sdmVycy9fX2luaXRfXy5weQ==..c7d429d47f2e44566ad76e30b8977572b074db44_Zmx1aWRzaW0vc29sdmVycy9fX2luaXRfXy5weQ== 100644
--- a/fluidsim/solvers/__init__.py
+++ b/fluidsim/solvers/__init__.py
@@ -18,4 +18,11 @@
    sphere
    models0d
 
+Provides:
+
+.. autosummary::
+   :toctree:
+
+   pkgload
+
 """
@@ -21,1 +28,19 @@
 """
+import numpy as _np
+
+
+def pkgload():
+    """Populate ``fluidsim.solvers`` package namespace with the solvers for easier
+    import.
+
+    """
+    from ..util.util import available_solver_keys
+    solvers = available_solver_keys()
+
+    solver_pkgs = (f"fluidsim.solvers.{solver}"
+                   for solver in solvers)
+    solver_modules = (f"fluidsim.solvers.{solver}.solver"
+                      for solver in solvers)
+
+    _np.pkgload(*solver_pkgs)
+    _np.pkgload(*solver_modules)