The simple answer is that you would use inverse iteration (subspace or with deflation). This is basically the power method (repeatedly multiplying the matrix by a vector and normalizing, singling out the eigenvector corresponding to the eigenvalue of largest magnitude) applied to $A^{-1}$. Since you desire the $k$ closest to the origin, you need to use some kind of block method or orthogonalization, otherwise the power iteration converges to the same eigenvector each time.
Practically, these simple algorithms are too primitive to be useful (they would take too long). More sophisticated variations of these include the (block) Lanczos or Arnoldi or Jacobi-Davidson methods. Be careful to note the restriction on the matrix properties (e.g. classical Lanczos only works on Hermitian matrices). Note that since you want the smallest eigenvalues, you need to be able to apply $A^{-1}x$ in these iterations. So each iteration requires a sparse matrix solve. Your 3D spatial grid structure might be amenable for multigrid techniques to speed up the solve, or if it's of moderate size you can do an (incomplete) LU or Cholesky decomposition. You can also use these as a preconditioner for the solver, which typically the eigenvalue methods take as a parameter.
Practically, Matlab has eigs, and this free book has a decent listing of existing code.
