We can obtain a solution by modifying the above PDE. To do so, we introduce a smoothing function sm[x] which effectively removes the problematic functions outside of the potential. Then, thanks to @xzczd 's suggestions and the links in the comments, we implement the absorbing boundary condition (abc) for the simple wave equation.
Geometry:
{xin, xmax, tmax, pts} = {-60, 100, 300, 350}; {x0, xp, w, wp, v} = {50, 0, 2.5, 2, -1}; c = 0.4; op = xmax - 1; g[x_, t_] := Exp[-(((x - x0) - v*t)/w)^2] gp[x_] := Exp[-((x - xp)/wp)^2] sm[x_] := (-Tanh[2 wp x] + 1)/2 Wave equation and boundary conditions:
mol[n_Integer] := {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n}} weqn = D[psi[x, t], {x, 2}] - D[psi[x, t], {t, 2}] - 2 c*I*D[psi[x, t], {t, 1}] sm[x] + (c^2 sm[x] - gp[x])* psi[x, t] == 0; ic = {psi[x, 0] == g[x, 0], Derivative[0, 1][psi][x, 0] == Derivative[0, 1][g][x, 0]}; bcd = {psi[xin, t] == 0}; abcr = D[psi[x, t], x] + D[psi[x, t], t] == 0 /. x -> xmax; Solution:
nsol = NDSolveValue[{weqn, ic, bcd, abcr}, psi, {x, xin, xmax}, {t, 0, tmax}, Method -> mol[pts]];