The thermodynamic properties of spin systems are evaluated with Monte Carlo methods. A review of classical thermodynamics is followed by a discussion of critical exponents. The Monte Carlo method is then applied to the two-dimensional Ising model with the goal of determining the phase diagram for magnetization. Boundary conditions, the reweighting method, autocorrelation, and critical slowing down are all explored. Cluster algorithms for overcoming critical slowing down are developed next and shown to dramatically reduce autocorrelation. A variety of spin systems that illustrate first, second, and infinite order (topological) phase transitions are explored. Finally, applications to random systems called spin glasses and to neural networks are briefly reviewed.