docs: def _initialize_stretching_matrix(self):
          """ Set up the stretching matrix """
  
          self.S = np.zeros((self.nz, self.nz))
  
          if (self.nz==2) and (self.rd) and (self.delta):
  
              self.del1 = self.delta/(self.delta+1.)
              self.del2 = (self.delta+1.)**-1
              self.Us = self.Ubg[0]-self.Ubg[1]
  
              self.F1 = self.rd**-2 / (1.+self.delta)
              self.F2 = self.delta*self.F1
              self.S[0,0], self.S[0,1] = -self.F1,  self.F1
              self.S[1,0], self.S[1,1] =  self.F2, -self.F2
  
          else:
  
              for i in range(self.nz):
  
                  if i == 0:
                      self.S[i,i]   = -self.f2/self.Hi[i]/self.gpi[i]
                      self.S[i,i+1] =  self.f2/self.Hi[i]/self.gpi[i]
  
                  elif i == self.nz-1:
                      self.S[i,i]   = -self.f2/self.Hi[i]/self.gpi[i-1]
                      self.S[i,i-1] =  self.f2/self.Hi[i]/self.gpi[i-1]
  
                  else:
                      self.S[i,i-1] = self.f2/self.Hi[i]/self.gpi[i-1]
                      self.S[i,i]   = -(self.f2/self.Hi[i]/self.gpi[i] +
                                          self.f2/self.Hi[i]/self.gpi[i-1])
                      self.S[i,i+1] = self.f2/self.Hi[i]/self.gpi[i]