In recent years, the characteristics of VCSELs have improved enormously that they have attracted more considerable interest. The corresponding developments of circuit-level VCSEL models are required for use in the design and simulation of optoelectronic applications. Unfortunately, existing models lack either the computational efficiency or the comprehensiveness warranted. In this paper, we present a VCSEL model that addresses the spatial dependent and thermal behavior of VCSELs based on multimode rate equations without sacrificing the numerical efficiency demanded by the circuit-level simulation of optoelectronic systems. Moreover, the thermal behavior is modeled by a laser diode based on Tucker's model with the parameters extracted directly by method of Gao and temperature dependent parasitic resistor and reverse saturation current are also considered. The equivalent circuit of the model is given and it is implemented into SPICE-like simulators to simulate the dc, ac and transient features of a 863 nm bottom-emitting AlGaAs VCSEL. The simulated results exhibit good agreements with references.