We describe a method for computing the evolution of a spherical stellar system using a fluid-dynamical approach based on the numerical solution of moment equations derived from the Boltzmann equation. Moments of the velocity distribution up to the fourth order are included in this treatment. In order to represent relaxation effects, `collision terms' are included and are evaluated using the Fokker-Planck equation. The principal assumption of the method is that the velocity distribution deviates by only a small amount from a Maxwellian distribution; however, it has been attempted to allow for a large anisotropy of the velocity distribution in the outer part of the system. To illustrate the application of the method, we apply it to the problem of the `gravothermal catastrophe' of an isothermal sphere, as recently discussed by Lynden-Bell & Wood (1968). The results confirm the predictions of LyndenBell & Wood concerning the instability and thermal runaway of an enclosed isothermal sphere when the central condensation becomes too high; however, the time scale for the ensuing' gravothermal catastrophe `is found to be much longer than the relaxation time. A discussion of the results shows that the rate of evolution is intimately related to the overall structure of the system; the more closely it approaches an isothermal structure with a Maxwellian velocity distribution, the slower is the evolution.