code: def equation(self):
          
          mat = np.zeros((self.nunknowns, self.model.neq))
          rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
          rhs[0:self.nlayers - 1] = 0.0
          rhs[self.nlayers - 1] = self.Qc
          ieq = 0
          for e in self.model.elementlist:
              if e.nunknowns > 0:
                  head = e.potinflayers(self.xc, self.yc, self.layers) / self.aq.Tcol[self.layers, :]
                  mat[0:self.nlayers - 1, ieq:ieq + e.nunknowns] = head[:-1] - head[1:]
                  if e == self:
                      for i in range(self.nlayers - 1):
                          mat[i, ieq + i] -= self.resfac[i]
                          mat[i, ieq + i + 1] += self.resfac[i + 1]
                      mat[self.nlayers - 1, ieq:ieq + self.nlayers] = 1.0
                  ieq += e.nunknowns
              else:
                  head = e.potentiallayers(self.xc, self.yc, self.layers) / self.aq.T[self.layers]
                  rhs[0:self.nlayers - 1] -= head[:-1] - head[1:]
          return mat, rhs
docs: Mix-in class that returns matrix rows for mscreen condition.
          Mscreen condition applied at each control point separately (so not like in ditch).
          Returns matrix part (nunknowns,neq)
          Returns rhs part (nunknowns)