The conventional rate equations for a Q-switched laser are augmented to explicitly include the effects of time-and level-dependent pumping, thermalization among the sublevels in the upper and lower multiplets, and multiplet relaxation in a homogeneously broadened four-level laser medium. To make the numerical computations more generally valid, we introduce a number of dimensionless variables. We show that the initial set of five coupled differential equations can be reduced to a simple set of two coupled equations for the inversion density and photon flux. Via numerical modeling, we have investigated the manner in which both thermalization and lower multiplet relaxation affect Nd:YAG laser characteristics such as output energy and temporal waveform. Our numerical results confirm earlier predictions that the Q-switched Nd:YAG laser output energy increases monotonically by a factor of 3.33 as one progresses from the assumption of slow to rapid thermalization and by an additional factor of 1.46 if one further assumes a terminal multiplet relaxation which is fast relative to the resonator photon decay time. We also find that the laser pulsewidth is substantially broadened when the resonator photon decay time is comparable to the thermalization, and to a lesser extent, the terminal multiplet relaxation times. Furthermore, when the Q-switched laser is operated near threshold and the photon decay time is significantly less than the thermalization time constant, the laser will produce a damped pulse train resulting from periodic replenishment of the population inversion through thermalization among the Stark sublevels, Our energy and pulsewidth calculations reproduce the measured values of two well-characterized diode-pumped Nd:YAG slab lasers quite well provided one assumes a thermalization time among nondegenerate Stark sublevels of 3-5 ns, This result contrasts with other experimental work which suggests subnanosecond terminal multiplet lifetimes. Possible reasons for the discrepancies are explored.