The paper presents a method for the recursive solution of models with a continuum of heterogeneous agents. Following Krusell and Smith (1998) and others, it is assumed that the wealth distribution in the economy can be represented, to a sufficient degree of accuracy, by a finite number of statistics. The method then uses a discretization of the state space on a finite grid. The method combines shape-preserving interpolation of the value function in endogenous variables with a simplicial linear interpolation in nonlinear transformations of other state variables. A shape-preserving interpolation scheme is developed that is particularly suitable for value functions that arise from capital accumulation problems. This interpolation scheme leads to an algorithm that - is stable - achieves a high level of accuracy with a relatively low computational effort. A computational advantage of the method is that it can be easily parallelized. The method is illustrated by an application to a standard model of heterogeneous agents that face uninsurable income risk, with endogenous labor supply and aggregate as well as idiosyncratic risk. The accuracy of the solution is investigated by a technique developed in Reiter (2000) which measures the optimizing agents' value loss which arises from using the numerical rather than the exact solution. While the main focus of the paper is on the stability and accuracy of the method, it also pays attention to the efficient implementation of the algorithm and investigates acceleration schemes to speed up the backward iterations. It is shown that a reasonably accurate solution can be obtained on a personal computer in a few minutes. References: Krusell, P. and Smith, A.A., Jr., Income and Wealth Heterogeneity in the Macroeconomy, JPE Vol. 106, 867-96, 1998 Reiter, M.: Estimating the Accuracy of Numerical Solutions to Dynamic Optimization Problems, Universitat Pompeu Fabra, 2000
