  We report the production and benchmarking of several refinements of the power
method that enable the computation of multiple extremal eigenpairs of very
large matrices. In these refinements we used an observation by Booth that has
made possible the calculation of up to the 10$^{th}$ eigenpair for simple test
problems simulating the transport of neutrons in the steady state of a nuclear
reactor. Here, we summarize our techniques and efforts to-date on determining
mainly just the two largest or two smallest eigenpairs. To illustrate the
effectiveness of the techniques, we determined the two extremal eigenpairs of a
cyclic matrix, the transfer matrix of the two-dimensional Ising model, and the
Hamiltonian matrix of the one-dimensional Hubbard model.
