code: def fd_solve(self):
      
      
      if self.Debug:
        try:
          # Will fail if scalar
          print("self.Te", self.Te.shape)
        except:
          pass
        print("self.qs", self.qs.shape)
        self.calc_max_flexural_wavelength()
        print("maxFlexuralWavelength_ncells: (x, y):", self.maxFlexuralWavelength_ncells_x, self.maxFlexuralWavelength_ncells_y)
      
      q0vector = self.qs.reshape(-1, order='C')
      if self.Solver == "iterative" or self.Solver == "Iterative":
        if self.Debug:
          print("Using generalized minimal residual method for iterative solution")
        if self.Verbose:
          print("Converging to a tolerance of", self.iterative_ConvergenceTolerance, "m between iterations")
        wvector = scipy.sparse.linalg.isolve.lgmres(self.coeff_matrix, q0vector)#, tol=1E-10)#,x0=woldvector)#,x0=wvector,tol=1E-15)    
        wvector = wvector[0] # Reach into tuple to get my array back
      else:
        if self.Solver == "direct" or self.Solver == "Direct":
          if self.Debug:
            print("Using direct solution with UMFpack")
        else:
          if self.Quiet == False:
            print("Solution type not understood:")
            print("Defaulting to direct solution with UMFpack")
        wvector = scipy.sparse.linalg.spsolve(self.coeff_matrix, q0vector, use_umfpack=True)
  
      # Reshape into grid
      self.w = -wvector.reshape(self.qs.shape)
      self.w_padded = self.w.copy()
docs: w = fd_solve()
      Sparse flexural response calculation.
      Can be performed by direct factorization with UMFpack (defuault) 
      or by an iterative minimum residual technique
      These are both the fastest of the standard Scipy builtin techniques in 
      their respective classes 
      Requires the coefficient matrix from "2D.coeff_matrix"