def fd_solve(self):
    """
    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"
    """
    
    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()