diff --git a/.flake8 b/.flake8
new file mode 100644
index 0000000..317f180
--- /dev/null
+++ b/.flake8
@@ -0,0 +1,4 @@
+[flake8]
+max-line-length = 88
+extend-ignore = E203
+exclude = tests
diff --git a/src/matadi/__about__.py b/src/matadi/__about__.py
index 970659c..86205cb 100644
--- a/src/matadi/__about__.py
+++ b/src/matadi/__about__.py
@@ -1 +1 @@
-__version__ = "0.1.16"
+__version__ = "0.1.17"
diff --git a/src/matadi/__init__.py b/src/matadi/__init__.py
index 2aa593b..acbaeff 100644
--- a/src/matadi/__init__.py
+++ b/src/matadi/__init__.py
@@ -1,10 +1,5 @@
-import casadi
-
-from .__about__ import __version__
-
-Variable = casadi.SX.sym
-
 from . import math, models
+from .__about__ import __version__
 from ._lab_compressible import LabCompressible
 from ._lab_compressible import LabCompressible as Lab
 from ._lab_incompressible import LabIncompressible
@@ -26,7 +21,30 @@
     TwoFieldVariation,
     TwoFieldVariationPlaneStrain,
 )
+from ._variable import Variable
 
 __all__ = [
     "__version__",
+    "math",
+    "models",
+    "LabCompressible",
+    "LabIncompressible",
+    "Lab",
+    "Function",
+    "FunctionScalar",
+    "FunctionTensor",
+    "Material",
+    "MaterialScalar",
+    "MaterialTensor",
+    "MaterialComposite",
+    "MaterialHyperelastic",
+    "MaterialHyperelasticPlaneStrain",
+    "MaterialHyperelasticPlaneStressIncompressible",
+    "MaterialHyperelasticPlaneStressLinearElastic",
+    "MaterialTensorGeneral",
+    "ThreeFieldVariation",
+    "ThreeFieldVariationPlaneStrain",
+    "TwoFieldVariation",
+    "TwoFieldVariationPlaneStrain",
+    "Variable",
 ]
diff --git a/src/matadi/_lab_incompressible.py b/src/matadi/_lab_incompressible.py
index 01922b9..4336515 100644
--- a/src/matadi/_lab_incompressible.py
+++ b/src/matadi/_lab_incompressible.py
@@ -2,7 +2,6 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy.optimize import root
 
 
 class LabIncompressible:
diff --git a/src/matadi/_material.py b/src/matadi/_material.py
index bb7cc76..735aa2f 100644
--- a/src/matadi/_material.py
+++ b/src/matadi/_material.py
@@ -3,8 +3,8 @@
 import casadi as ca
 import numpy as np
 
-from . import Variable
 from ._apply import apply
+from ._variable import Variable
 
 
 class Function:
diff --git a/src/matadi/_templates.py b/src/matadi/_templates.py
index 5291022..2622f2c 100644
--- a/src/matadi/_templates.py
+++ b/src/matadi/_templates.py
@@ -1,7 +1,7 @@
 import numpy as np
 
-from . import Variable
 from ._material import Material, MaterialTensor
+from ._variable import Variable
 from .math import Function, det, eye
 from .math import gradient as grad
 from .math import horzcat, trace, vertcat, zeros
@@ -241,7 +241,7 @@ def __init__(self, fun, statevars_shape=(1, 1), x=None, triu=True, **kwargs):
         try:
             # displacement-pressure split
             x.append(fun.p)
-        except:
+        except AttributeError:
             pass
 
         # add state variables
diff --git a/src/matadi/_variable.py b/src/matadi/_variable.py
new file mode 100644
index 0000000..6fd43a7
--- /dev/null
+++ b/src/matadi/_variable.py
@@ -0,0 +1,3 @@
+import casadi
+
+Variable = casadi.SX.sym
diff --git a/src/matadi/math.py b/src/matadi/math.py
index f68e5a8..e952d3c 100644
--- a/src/matadi/math.py
+++ b/src/matadi/math.py
@@ -67,6 +67,73 @@
     vertsplit,
 )
 
+__all__ = [
+    "DM",
+    "MX",
+    "SX",
+    "Function",
+    "acos",
+    "acosh",
+    "adj",
+    "asin",
+    "asinh",
+    "atan",
+    "atan2",
+    "atanh",
+    "ceil",
+    "cofactor",
+    "cos",
+    "cosh",
+    "cross",
+    "det",
+    "diag",
+    "_dot",
+    "eig_symbolic",
+    "erf",
+    "erfinv",
+    "exp",
+    "fabs",
+    "find",
+    "floor",
+    "fmax",
+    "fmin",
+    "gradient",
+    "hessian",
+    "horzcat",
+    "horzsplit",
+    "if_else",
+    "inv",
+    "ldl",
+    "linspace",
+    "log",
+    "logic_and",
+    "logic_not",
+    "logic_or",
+    "mmax",
+    "mmin",
+    "norm_1",
+    "pi",
+    "qr",
+    "repmat",
+    "reshape",
+    "sign",
+    "sin",
+    "sinh",
+    "sqrt",
+    "sum1",
+    "sum2",
+    "sumsqr",
+    "tan",
+    "tanh",
+    "times",
+    "trace",
+    "transpose",
+    "tril",
+    "triu",
+    "vertcat",
+    "vertsplit",
+]
+
 eye = SX.eye
 ones = SX.ones
 zeros = SX.zeros
@@ -134,11 +201,11 @@ def tresca(C):
 def mexp(C, eps=8e-5):
     "Exponential Function of a Matrix."
     w = eigvals(C, eps=eps)
-    I = SX.eye(3)
+    eye = SX.eye(3)
 
-    M1 = (C - w[1] * I) * (C - w[2] * I) / (w[0] - w[1]) / (w[0] - w[2])
-    M2 = (C - w[2] * I) * (C - w[0] * I) / (w[1] - w[2]) / (w[1] - w[0])
-    M3 = (C - w[0] * I) * (C - w[1] * I) / (w[2] - w[0]) / (w[2] - w[1])
+    M1 = (C - w[1] * eye) * (C - w[2] * eye) / (w[0] - w[1]) / (w[0] - w[2])
+    M2 = (C - w[2] * eye) * (C - w[0] * eye) / (w[1] - w[2]) / (w[1] - w[0])
+    M3 = (C - w[0] * eye) * (C - w[1] * eye) / (w[2] - w[0]) / (w[2] - w[1])
 
     return exp(w[0]) * M1 + exp(w[1]) * M2 + exp(w[2]) * M3
 
diff --git a/src/matadi/models/__init__.py b/src/matadi/models/__init__.py
index 4840a8e..f3d6d5a 100644
--- a/src/matadi/models/__init__.py
+++ b/src/matadi/models/__init__.py
@@ -22,3 +22,30 @@
 from ._templates import Morph, NeoHookeOgdenRoxburgh, Viscoelastic
 from ._viscoelasticity import finite_strain_viscoelastic
 from .microsphere.nonaffine import miehe_goektepe_lulei
+
+__all__ = [
+    "microsphere",
+    "displacement_pressure_split",
+    "isochoric_volumetric_split",
+    "volumetric",
+    "fiber",
+    "fiber_family",
+    "holzapfel_gasser_ogden",
+    "arruda_boyce",
+    "extended_tube",
+    "linear_elastic",
+    "mooney_rivlin",
+    "neo_hooke",
+    "ogden",
+    "saint_venant_kirchhoff",
+    "third_order_deformation",
+    "van_der_waals",
+    "yeoh",
+    "morph",
+    "ogden_roxburgh",
+    "Morph",
+    "NeoHookeOgdenRoxburgh",
+    "Viscoelastic",
+    "finite_strain_viscoelastic",
+    "miehe_goektepe_lulei",
+]
diff --git a/src/matadi/models/_helpers.py b/src/matadi/models/_helpers.py
index 1e03be6..80f79a6 100644
--- a/src/matadi/models/_helpers.py
+++ b/src/matadi/models/_helpers.py
@@ -1,7 +1,7 @@
 from copy import deepcopy
 from functools import wraps
 
-from .. import Variable
+from .._variable import Variable
 from ..math import cof, det, gradient, trace
 
 
diff --git a/src/matadi/models/_hyperelasticity_anisotropic.py b/src/matadi/models/_hyperelasticity_anisotropic.py
index 4255086..4c7d7c7 100644
--- a/src/matadi/models/_hyperelasticity_anisotropic.py
+++ b/src/matadi/models/_hyperelasticity_anisotropic.py
@@ -1,7 +1,6 @@
 from ..math import (
     DM,
     cos,
-    det,
     exp,
     if_else,
     invariants,
diff --git a/src/matadi/models/_hyperelasticity_isotropic.py b/src/matadi/models/_hyperelasticity_isotropic.py
index 9278772..1ca1258 100644
--- a/src/matadi/models/_hyperelasticity_isotropic.py
+++ b/src/matadi/models/_hyperelasticity_isotropic.py
@@ -22,7 +22,7 @@ def neo_hooke(F, C10):
 
 
 @isochoric_volumetric_split
-def mooney_rivlin(F, C10, C01):
+def mooney_rivlin(F, C10, C01=0):
     C = transpose(F) @ F
     I1 = trace(C)
     I2 = (trace(C) ** 2 - trace(C @ C)) / 2
@@ -30,7 +30,7 @@ def mooney_rivlin(F, C10, C01):
 
 
 @isochoric_volumetric_split
-def yeoh(F, C10, C20, C30):
+def yeoh(F, C10, C20=0, C30=0):
     J = det(F)
     C = transpose(F) @ F
     I1 = J ** (-2 / 3) * trace(C)
@@ -38,7 +38,7 @@ def yeoh(F, C10, C20, C30):
 
 
 @isochoric_volumetric_split
-def third_order_deformation(F, C10, C01, C11, C20, C30):
+def third_order_deformation(F, C10, C01=0, C11=0, C20=0, C30=0):
     C = transpose(F) @ F
     I1 = trace(C)
     I2 = (trace(C) ** 2 - trace(C @ C)) / 2
@@ -96,8 +96,8 @@ def van_der_waals(F, mu, limit, a, beta):
     C = transpose(F) @ F
     I1 = trace(C)
     I2 = (trace(C) ** 2 - trace(C @ C)) / 2
-    I = (1 - beta) * I1 + beta * I2
-    eta = sqrt((I - 3) / (limit**2 - 3))
+    Im = (1 - beta) * I1 + beta * I2
+    eta = sqrt((Im - 3) / (limit**2 - 3))
     return mu * (
-        -(limit**2 - 3) * (log(1 - eta) + eta) - 2 / 3 * a * ((I - 3) / 2) ** (3 / 2)
+        -(limit**2 - 3) * (log(1 - eta) + eta) - 2 / 3 * a * ((Im - 3) / 2) ** (3 / 2)
     )
diff --git a/src/matadi/models/_misc.py b/src/matadi/models/_misc.py
index abf5e32..f194e5f 100644
--- a/src/matadi/models/_misc.py
+++ b/src/matadi/models/_misc.py
@@ -51,7 +51,9 @@ def morph(x, p1, p2, p3, p4, p5, p6, p7, p8):
     CT_CTS = if_else(CTS > 0, CT / CTS, CT)
 
     # MORPH deformation-dependent material parameters
-    f = lambda x: 1 / sqrt(1 + x**2)
+    def f(x):
+        return 1 / sqrt(1 + x**2)
+
     a = p1 + p2 * f(p3 * CTS)
     b = p4 * f(p3 * CTS)
     c = p5 * CTS * (1 - f(CTS / p6))
diff --git a/src/matadi/models/microsphere/__init__.py b/src/matadi/models/microsphere/__init__.py
index fee0416..79ae8b6 100644
--- a/src/matadi/models/microsphere/__init__.py
+++ b/src/matadi/models/microsphere/__init__.py
@@ -1,2 +1,12 @@
 from . import affine, nonaffine, quadrature
 from ._chain import gauss, langevin, langevin2, linear
+
+__all__ = [
+    "affine",
+    "nonaffine",
+    "quadrature",
+    "gauss",
+    "langevin",
+    "langevin2",
+    "linear",
+]
diff --git a/src/matadi/models/microsphere/affine/__init__.py b/src/matadi/models/microsphere/affine/__init__.py
index 9c2cb47..39f3fef 100644
--- a/src/matadi/models/microsphere/affine/__init__.py
+++ b/src/matadi/models/microsphere/affine/__init__.py
@@ -1,3 +1,9 @@
 from ._models import microsphere_affine_force as force
 from ._models import microsphere_affine_stretch as stretch
 from ._models import microsphere_affine_tube as tube
+
+__all__ = [
+    "force",
+    "stretch",
+    "tube",
+]
diff --git a/src/matadi/models/microsphere/nonaffine/__init__.py b/src/matadi/models/microsphere/nonaffine/__init__.py
index 9a1d577..6caa501 100644
--- a/src/matadi/models/microsphere/nonaffine/__init__.py
+++ b/src/matadi/models/microsphere/nonaffine/__init__.py
@@ -1,3 +1,9 @@
 from ._models import microsphere_nonaffine_miehe_goektepe_lulei as miehe_goektepe_lulei
 from ._models import microsphere_nonaffine_stretch as stretch
 from ._models import microsphere_nonaffine_tube as tube
+
+__all__ = [
+    "miehe_goektepe_lulei",
+    "stretch",
+    "tube",
+]
diff --git a/src/matadi/models/microsphere/quadrature/__init__.py b/src/matadi/models/microsphere/quadrature/__init__.py
index 597df36..c604042 100644
--- a/src/matadi/models/microsphere/quadrature/__init__.py
+++ b/src/matadi/models/microsphere/quadrature/__init__.py
@@ -1 +1,5 @@
 from ._bazant_oh import BazantOh
+
+__all__ = [
+    "BazantOh",
+]