In this paper we solve the benchmark heterogeneous agents model by Aiyagari (1994) using Monte Carlo methods. In addition, the idiosyncratic shocks process is approximated using Tauchen's (1986) method. This we go beyond the 2 by 2 Markov matrix approximation of the AR(1) stochastic process. The code is written in MATLAB. The computation time is much faster than the one written by Heer and Maussner (2008) in FORTRAN. This model also solves Mehra-Prescott's puzzle and generates a risk-free interest rate that is much closer to the one we observe in data.