Research on the development of numerical techniques to simulate the space-time evolution of large tokamak plasmas is reported. A nonuniform spatial mesh technique is employed to allow more accurate calculations in the boundary of reactor-size plasmas. A box integration method is used to maintain the accuracy of central differencing on the nonuniform spatial mesh and to preserve both the particle and energy flux. A variable implicit technique is used for the time expansion. The time-centered (Crank-Nicholson) technique used in most other models generally offers greater accuracy but can lead to severe limitations on the time step. Somewhat more implicit treatments can remove the numerical limitations on the time step without seriously affecting accuracy. The physical time scales, which can change by several orders of magnitude from startup to equilibrium, can then be used to continually adjust the time step throughout a calculation. Sample calculations are presented for a near-term tokamak engineering test reactor and a conceptual tokamak power reactor, UWMAK-III.