code: def fd_solve(self):
      
      
      if self.Debug:
        print("qs", self.qs.shape)
        print("Te", self.Te.shape)
        self.calc_max_flexural_wavelength()
        print("maxFlexuralWavelength_ncells', self.maxFlexuralWavelength_ncells")
      
      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")
        # qs negative so bends down with positive load, bends up with neative load 
        # (i.e. material removed)
        w = isolve.lgmres(self.coeff_matrix, -self.qs, tol=self.iterative_ConvergenceTolerance)  
        self.w = w[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:
          print("Solution type not understood:")
          print("Defaulting to direct solution with UMFpack")
        # UMFpack is now the default, but setting true just to be sure in case
        # anything changes
        # qs negative so bends down with positive load, bends up with neative load 
        # (i.e. material removed)
        self.w = spsolve(self.coeff_matrix, -self.qs, use_umfpack=True)
      
      if self.Debug:
        print("w.shape:")
        print(self.w.shape)
        print("w:")
        print(self.w)
docs: w = fd_solve()
        where coeff is the sparse coefficient matrix output from function
        coeff_matrix and qs is the array of loads (stresses)
  
      Sparse solver for one-dimensional flexure of an elastic plate