docs: def solve(self, neigs=4, tol=0, guess=None, mode_profiles=True, initial_mode_guess=None):
          """
          This function finds the eigenmodes.
  
          Parameters
          ----------
          neigs : int
              number of eigenmodes to find
          tol : float
              Relative accuracy for eigenvalues. The default value of 0 implies machine precision.
          guess : float
              a guess for the refractive index. Only finds eigenvectors with an effective refractive index
              higher than this value.
  
          Returns
          -------
          self : an instance of the VFDModeSolver class
              obtain the fields of interest for specific modes using, for example:
              solver = EMpy.modesolvers.FD.VFDModeSolver(wavelength, x, y, epsf, boundary).solve()
              Ex = solver.modes[0].Ex
              Ey = solver.modes[0].Ey
              Ez = solver.modes[0].Ez
          """
  
          from scipy.sparse.linalg import eigen
  
          self.nmodes = neigs
          self.tol = tol
  
          A = self.build_matrix()
  
          if guess is not None:
              # calculate shift for eigs function
              k = 2 * numpy.pi / self.wl
              shift = (guess * k) ** 2
          else:
              shift = None
  
          [eigvals, eigvecs] = eigen.eigs(A,
                                          k=neigs,
                                          which='LR',
                                          tol=0.001,
                                          ncv=None,
                                          v0 = initial_mode_guess,
                                          return_eigenvectors=mode_profiles,
                                          sigma=shift)
  
          neffs = self.wl * scipy.sqrt(eigvals) / (2 * numpy.pi)
          if mode_profiles:
              Hxs = []
              Hys = []
              nx = self.nx
              ny = self.ny
              for ieig in range(neigs):
                  Hxs.append(eigvecs[:nx * ny, ieig].reshape(nx, ny))
                  Hys.append(eigvecs[nx * ny:, ieig].reshape(nx, ny))
  
          # sort the modes
          idx = numpy.flipud(numpy.argsort(neffs))
          neffs = neffs[idx]
          self.neff = neffs
          if mode_profiles:
              tmpx = []
              tmpy = []
              for i in idx:
                  tmpx.append(Hxs[i])
                  tmpy.append(Hys[i])
              Hxs = tmpx
              Hys = tmpy
  
              [Hzs, Exs, Eys, Ezs] = self.compute_other_fields(neffs, Hxs, Hys)
  
              self.modes = []
              for (neff, Hx, Hy, Hz, Ex, Ey, Ez) in zip(neffs, Hxs, Hys, Hzs, Exs, Eys, Ezs):
                  self.modes.append(
                      FDMode(self.wl, self.x, self.y, neff, Ey, Ex, Ez, Hy, Hx, Hz).normalize())
  
          return self