  We introduce the Chebyshev pseudosite matrix product state approach
(ChePSMPS) as a solver for cluster perturbation theory (CPT), crucial for
simulating spectral functions in two-dimensional electron-phonon ($e$-ph)
coupling systems. ChePSMPS distinguishes itself from conventional exact
diagonalization solvers by supporting larger clusters, thereby significantly
mitigating finite-size effects. Free from the fermion sign problem, ChePSMPS
enhances its ability to explore $e$-ph effects and generate high-resolution
spectral functions in doped Mott insulators. We use this method to simulate the
spectra for both one- and two-dimensional Hubbard-Holstein models, highlighting
its superiority over other methods. Our findings validate ChePSMPS as a
powerful and reliable Green's function solver. In conjunction with embedding
methods, ChePSMPS emerges as an essential tool for simulating strongly
correlated $e$-ph coupling systems.
