docs: def equation(self):
          """Mix-in class that returns matrix rows for difference in head between inside and
          outside equals zeros
          Returns matrix part (nunknowns,neq)
          Returns rhs part nunknowns
          """
          mat = np.empty((self.nunknowns, self.model.neq))
          rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
          for icp in range(self.ncp):
              istart = icp * self.nlayers
              ieq = 0
              for e in self.model.elementlist:
                  if e.nunknowns > 0:
                      qxin, qyin = e.disvecinflayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin)
                      qxout, qyout = e.disvecinflayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout)
                      mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                          (qxin - qxout) * self.cosnorm[icp] + (qyin - qyout) * self.sinnorm[icp]
                      ieq += e.nunknowns
                  else:
                      qxin, qyin = e.disveclayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin)
                      qxout, qyout = e.disveclayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout)
                      rhs[istart:istart + self.nlayers] -= (qxin - qxout) * self.cosnorm[icp] + (qyin - qyout) * \
                                                           self.sinnorm[icp]
          return mat, rhs