Shale is a common rock in oil and gas extraction, and the study of its nonlinear mechanical behavior is crucial for the development of engineering techniques such as hydraulic fracturing. This paper establishes a new coupled elastic–plastic damage model based on the second law of thermodynamics, the strain equivalence principle, the non-associated flow rule, and the Drucker–Prager yield criterion. This model is used to describe the mechanical behavior of shale before and after peak strength and has been implemented in ABAQUS via UMAT for numerical computation. The model comprehensively considers the quasi-brittle and anisotropic characteristics of shale, as well as the strength degradation caused by damage during both the elastic and plastic phases. A damage yield function has been established as a criterion for damage occurrence, and the constitutive integration algorithm has been derived using a regression mapping algorithm. Compared with experimental data from La Biche shale in Canada, the theoretical model accurately simulated the stress–strain curves and volumetric–axial strain curves of shale under confining pressures of 5 MPa, 25 MPa, and 50 MPa. When compared with experimental data from shale in Western Hubei and Eastern Chongqing, China, the model precisely fitted the stress–strain curves of shale at pressures of 30 MPa, 50 MPa, and 70 MPa, and at bedding angles of 0°, 22.5°, 45°, and 90°. This proves that the model can effectively predict the failure behavior of shale under different confining pressures and bedding angles. Additionally, a sensitivity analysis has been performed on parameters such as the plastic hardening rate b, damage evolution rate B_ω, weighting factor r, and damage softening parameter a. This research is expected to provide theoretical support for the efficient extraction technologies of shale oil and gas.