Via adopting the matrix exponential (ME) approach, efficient and tight complex-frequency-shifted perfectly matched layer (CFS-PML) formulations using the matched Z-transform (MZT) method are well organized and proposed to truncate the monolayer black phosphorous (BP) metasurface problems in the infrared regime. Due to introducing ME approach, the proposed MZT-ME-PML develops into the compact first-order differential matrix form used to the circumvent convolution manipulations, the variable replacements, or the formulation rearrangements during the mathematical derivations. Furthermore, the MZT-ME-PML can not only flexibly truncate arbitrary materials due to using the DB constitutive relation, but efficiently attenuate evanescent waves and reduce late-time reflections with the CFS scheme. Besides, in view of the discretization effect on Terahertz (THz) BP metasurface implementation, the MZT technique is added deliberately into this proposal in that it is more suitable for simulating fine structures. Finally, the ME-based proposal can maintain high absorption when it comes closer to terminate finite-difference time-domain (FDTD) domains so that smaller physical regions can be applied, leading to palpable improvements in memory and CPU time. The three-dimensional (3D) numerical simulations on monolayer BP THz problems have been carried out to illustrate the validity and accuracy of the proposal.