diff --git a/docs/03examples/wells_in_different_systems.ipynb b/docs/03examples/wells_in_different_systems.ipynb
index e6e6b4c..ac0c871 100644
--- a/docs/03examples/wells_in_different_systems.ipynb
+++ b/docs/03examples/wells_in_different_systems.ipynb
@@ -237,6 +237,154 @@
     "plt.title(\"note that heads are not computed at the same times\")\n",
     "plt.grid()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Well with entry resistance"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "self.neq  3\n",
+      "solution complete\n"
+     ]
+    }
+   ],
+   "source": [
+    "k = 25  # m/d\n",
+    "H = 20  # m\n",
+    "Ss = 1e-4 / H  # m^(-1)\n",
+    "rw = 0.3  # m\n",
+    "res = 0.1  # d\n",
+    "Q = 500  # m^3/d\n",
+    "ml = ttim.Model3D(kaq=k, z=np.linspace(H, 0, 6), Saq=Ss, tmin=0.001, tmax=10)\n",
+    "w = ttim.Well(ml, tsandQ=[(0, Q)], rw=rw, res=res, layers=[0, 1, 2])\n",
+    "ml.solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "head inside well at t=2:      [-4.18893381 -4.18893378 -4.18893388]\n",
+      "head just outside well at t=2 [-2.01308246 -1.99455564 -1.92770745]\n"
+     ]
+    }
+   ],
+   "source": [
+    "hin = w.headinside(t=2)[:, 0]\n",
+    "hout = ml.head(rw, 0, t=2, layers=[0, 1, 2])[:, 0]\n",
+    "print(\"head inside well at t=2:     \", hin)\n",
+    "print(\"head just outside well at t=2\", hout)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "discharge of screens at t=2 [164.05532741 165.45221426 170.49245612]\n",
+      "discharge from bc at t=2:   [164.05532741 165.45221426 170.49245612]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"discharge of screens at t=2\", w.discharge(2)[:, 0])\n",
+    "Qcheck = 2 * np.pi * rw * ml.aq.Haq[:3] * (hout - hin) / res\n",
+    "print(\"discharge from bc at t=2:  \", Qcheck)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Well with different entry resistance per layer"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "self.neq  3\n",
+      "solution complete\n"
+     ]
+    }
+   ],
+   "source": [
+    "k = 25  # m/d\n",
+    "H = 20  # m\n",
+    "Ss = 1e-4 / H  # m^(-1)\n",
+    "rw = 0.3  # m\n",
+    "res = [0.05, 0.1, 0.15]  # d\n",
+    "Q = 500  # m^3/d\n",
+    "ml = ttim.Model3D(kaq=k, z=np.linspace(H, 0, 6), Saq=Ss, tmin=0.001, tmax=10)\n",
+    "w = ttim.Well(ml, tsandQ=[(0, Q)], rw=rw, res=res, layers=[0, 1, 2])\n",
+    "ml.solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "head inside well at t=2:      [-3.95216559 -3.9521661  -3.95216657]\n",
+      "head just outside well at t=2 [-2.42278217 -1.93402805 -1.62034002]\n"
+     ]
+    }
+   ],
+   "source": [
+    "hin = w.headinside(t=2)[:, 0]\n",
+    "hout = ml.head(rw, 0, t=2, layers=[0, 1, 2])[:, 0]\n",
+    "print(\"head inside well at t=2:     \", hin)\n",
+    "print(\"head just outside well at t=2\", hout)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "discharge of screens at t=2 [230.62558655 152.1640244  117.21038657]\n",
+      "discharge from bc at t=2:   [230.62558655 152.1640244  117.21038657]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"discharge of screens at t=2\", w.discharge(2)[:, 0])\n",
+    "Qcheck = 2 * np.pi * rw * ml.aq.Haq[:3] * (hout - hin) / res\n",
+    "print(\"discharge from bc at t=2:  \", Qcheck)"
+   ]
   }
  ],
  "metadata": {
diff --git a/ttim/equation.py b/ttim/equation.py
index 7fc43dd..0b527dd 100644
--- a/ttim/equation.py
+++ b/ttim/equation.py
@@ -74,16 +74,7 @@ def equation(self):
                     np.pi * self.rc**2 * self.model.p * head[0, :]
                 )
                 if e == self:
-                    disterm = (
-                        self.dischargeinflayers
-                        * self.res
-                        / (
-                            2
-                            * np.pi
-                            * self.rw
-                            * self.aq.Haq[self.layers][:, np.newaxis]
-                        )
-                    )
+                    disterm = self.dischargeinflayers * self.resfach[:, np.newaxis]
                     if self.nunknowns > 1:  # Multiple layers
                         for i in range(self.nunknowns - 1):
                             mat[i, ieq + i, :] -= disterm[i]
diff --git a/ttim/size_of_arrays.txt b/ttim/size_of_arrays.txt
index 20beefc..3334af6 100644
--- a/ttim/size_of_arrays.txt
+++ b/ttim/size_of_arrays.txt
@@ -5,4 +5,4 @@ self.potinf = self.nparam, aq.naq, self.model.npval
 self.potential = np.sum(self.parameters[:, :, np.newaxis, :] * self.potinf(x, y, aq), 1)
 ngvbc, naq, npval
 
-self.lab2.shape = (self.naq, self.model.nint, self.model.npint)
\ No newline at end of file
+self.lab2.shape = (self.naq, self.model.nint, self.model.npint)