We analyze numerical methods for the system coupling the (frequency averaged) radiation energy density Er to the matter energy balance equation. Since second order temporal differencing may give unphysical results, attention is focused on 1st order backward Euler temporal differencing and a Picard- Newton iterative scheme solves the nonlinear equations. We show that iterating on Er and the emission source B is more robust. Pseudo transient continuation (PsiTC) is introduced in order obtain a scheme that converges for wide ranges of initial conditions and time steps. We derive conditions for the PsiTC parameter guaranteeing physically meaningful, i.e., positive, (Er,B) iterates which can be bounded and shown to equilibrate. The scheme has been incorporated into a multiply- dimensioned, massively parallel, Eulerian, radiation-hydrodynamic code with AMR. We present three simulations: (1) Pomraning problem. (2) An ionization front propagating into tenuous hydrogen gas with a (nonlinear) Saha model for the equation of state. (3) Interacting, radiatively driven shocks propagating through real materials. **This work was performed under the auspices of the U.S. Dept. of Energy by the Univ. of Calif. Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.