Melting in a finite two-dimensional dusty plasma has recently been characterized experimentally [T. E. Sheridan, Phys. Plasmas 15, 103702 (2008)] in a system with ≈3900 dust particles. We model this experiment using the Metropolis algorithm to generate thermodynamic configurations for a two-dimensional system of identical particles confined in a parabolic well and interacting through a Debye (Yukawa) potential. Results are computed for a Debye shielding parameter a/λD=0.24, where a is the lattice constant and λD is the Debye length. Configurations are characterized using the pair and bond-orientational correlation functions, defect fractions, and correlation lengths. Distinct crystalline, hexatic, and liquid phases are observed. In the hexatic phase, the decay of the bond-orientational correlation goes as r−η6 where η6≈1, which is greater than the value η6=0.25 predicted by the theory of Kosterlitz, Thouless, Halperin, Nelson, and Young.