Large scale simulation of the colossal magnetoresistance effect in manganites using spin-fermion models is often hampered due to the high computational cost associated with computing the eigenvalues of each of the successive Hamiltonian matrices. Consequently, current spin-fermion model simulations contain no more than 6³ sites or the equivalent for lower dimensions. This imposes severe limitations on the kinds of physical systems that can be studied; for example, the Mn spin concentration in diluted semiconductors has to be high enough to be numerically tractable, and the study of many band systems becomes computationally difficult. This study presents an algorithm that directly updates the spectrum of a successive Hamiltonian matrix on the basis of the spectrum of the previous Hamiltonian matrix. This eigenvalue updating algorithm significantly reduces the computational bottleneck involved in recomputing the spectrum of the Hamiltonian matrices each time a local configurational change is accepted, thereby allowing the simulation of much larger lattice system sizes. The serial version of the algorithm is an order of magnitude faster than the approaches based on direct diagonalization. In addition, this algorithm is amenable to parallel computation and retains excellent accuracy even after many updates.