docs: def equation(self):
          """Mix-in class that returns matrix rows for potential-specified conditions.
          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
          rhs = self.pc.copy()
          for icp in range(self.ncp):
              istart = icp * self.nlayers
              # rhs[istart:istart+self.nlayers] = self.pc[]
              ieq = 0
              for e in self.model.elementlist:
                  if e.nunknowns > 0:
                      mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                          e.potinflayers(self.xc[icp], self.yc[icp], self.layers)
                      if e == self:
                          mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] -= self.resfac[icp]
                      ieq += e.nunknowns
                  else:
                      rhs[istart:istart + self.nlayers] -= \
                          e.potentiallayers(self.xc[icp], self.yc[icp],
                                            self.layers)  # Pretty cool that this works, really
          return mat, rhs