diff --git a/timflow/bessel/besselnumba.py b/timflow/bessel/besselnumba.py
index 2d9c13d..3d648b9 100644
--- a/timflow/bessel/besselnumba.py
+++ b/timflow/bessel/besselnumba.py
@@ -22,6 +22,7 @@
 
 fac = 1.0
 
+# max number of terms fixed to 20
 nrange = np.arange(21, dtype=np.float64)
 
 a = np.zeros(21, dtype=np.float64)
diff --git a/timflow/steady/inhomogeneity.py b/timflow/steady/inhomogeneity.py
index 7744d4c..a9021cc 100644
--- a/timflow/steady/inhomogeneity.py
+++ b/timflow/steady/inhomogeneity.py
@@ -424,7 +424,7 @@ def __repr__(self):
             "BuildingPit: layers "
             + str(list(self.layers))
             + ", "
-            + str(list(self.x, self.y))
+            + str(list(zip(self.x.tolist(), self.y.tolist(), strict=False)))
         )
 
     def isinside(self, x, y):
@@ -687,7 +687,9 @@ def __init__(
             layers = [0]
         if z is None:
             z = [1, 0]
-        (kaq, c, npor, ltype) = param_3d(kaq, z, kzoverkh, npor, topboundary, topres)
+        kaq, kzoverkh, c, npor, ltype = param_3d(
+            kaq, z, kzoverkh, npor, topboundary, topres
+        )
         if topboundary == "semi":
             z = np.hstack((z[0] + topthick, z))
         super().__init__(
@@ -703,6 +705,7 @@ def __init__(
             ndeg=ndeg,
             layers=layers,
         )
+        self.kzoverkh = kzoverkh  # add kzoverkh to inhomogeneity object
 
 
 class LeakyBuildingPit(BuildingPit):
@@ -812,7 +815,7 @@ def __repr__(self):
             "LeakyBuildingPit: layers "
             + str(list(self.layers))
             + ", "
-            + str(list(self.x, self.y))
+            + str(list(zip(self.x.tolist(), self.y.tolist(), strict=False)))
         )
 
     def create_elements(self):
@@ -1060,7 +1063,9 @@ def __init__(
             layers = [0]
         if z is None:
             z = [1, 0]
-        (kaq, c, npor, ltype) = param_3d(kaq, z, kzoverkh, npor, topboundary, topres)
+        kaq, kzoverkh, c, npor, ltype = param_3d(
+            kaq, z, kzoverkh, npor, topboundary, topres
+        )
         if topboundary == "semi":
             z = np.hstack((z[0] + topthick, z))
         super().__init__(
@@ -1077,6 +1082,7 @@ def __init__(
             layers=layers,
             res=res,
         )
+        self.kzoverkh = kzoverkh  # add kzoverkh to inhomogeneity object
 
 
 class AreaSinkInhom(Element):