Context: This question is relevant to the physical problem of an excited leaky cavity. The boundary condition on the inside of the cavity is Dirichlet. The other end of the cavity is partially blocked by a potential, and beyond that waves are free to escape to infinity. The dynamics are governed by a modified wave equation.

First set the geometry of the problem, with an example Gaussian perturbation and potential:

{xin, xmax, tmax} = {-300, 400, 3000};
{x0, xp, w, wp, v} = {-150, 0, 2, 8, 1};
c = 0.4;
g[x_,t_] := Exp[-(((x - x0) - v*t)/w)^2]
gp[x_] := Exp[-((x - xp)/wp)^2]

And define the modified wave equation:

weqn = D[psi[x, t], {x, 2}] - D[psi[x, t], {t, 2}] - 
    2 c*I*D[psi[x, t], {t, 1}] + (c^2 - gp[x])*
     psi[x, t] == -NeumannValue[-Derivative[0, 1][psi][x, t], 
     x == xmax];

Because waves should be purely outgoing outside of the potential, to avoid integrating out to large x we can implement an absorbing boundary condition by putting in the NeumannValue[...] term above for some xmax sufficiently far from the potential.

Initial conditions:

ic = {psi[x, 0] == g[x, 0], 
     Derivative[0, 1][psi][x, 0] == Derivative[0, 1][g][x, 0]};

Boundary condition for the inside of the cavity:

bcd = {psi[xin, t] == 0};

Attempt at a solution:

The absorbing boundary condition seems to require a FEM method, so we write

sol = psi /. 
   NDSolve[{weqn, ic, bcd}, psi, {x, xin, xmax}, {t, 0, tmax}, 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "InterpolationOrder" -> {psi -> 2}, 
         "MeshOptions" -> {"MaxCellMeasure" -> 1}}}][[1, 1]];

This works when c = 0, for which weqn is a simple wave equation with a potential. But when c = 0.4 the following warnings are generated,

NDSolve::femdpop: The NDSolve`FEM`FEMStiffnessElements operator failed.

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

The problem seems to be due to the third (imaginary) term in the PDE. So a natural fix would be to split the equation into its real and imaginary parts and solve the coupled PDE, so that all terms are real. In other words, we take psi = psiR + I*psiI and thus obtain

weqn1 = D[psiR[x, t], {x, 2}] - D[psiR[x, t], {t, 2}] + 
    2*c*D[psiI[x, t], {t, 1}] + (c^2 - gp[x])*
     psiR[x, t] == -NeumannValue[-Derivative[0, 1][psiR][x, t], 
     x == xmax];
weqn2 = D[psiI[x, t], {x, 2}] - D[psiI[x, t], {t, 2}] - 
    2*c*D[psiR[x, t], {t, 1}] + (c^2 - gp[x])*
     psiI[x, t] == -NeumannValue[-Derivative[0, 1][psiI][x, t], 
     x == xmax];
ic1 = {psiR[x, 0] == g[x, 0], 
   Derivative[0, 1][psiR][x, 0] == Derivative[0, 1][g][x, 0]};
ic2 = {psiI[x, 0] == 0, Derivative[0, 1][psiI][x, 0] == 0};
bcd1 = {psiR[xin, t] == 0};
bcd2 = {psiI[xin, t] == 0};

Now we solve

solP = NDSolve[{weqn1, weqn2, ic1, ic2, bcd1, bcd2}, {psiR, psiI}, {x,
      xin, xmax}, {t, 0, tmax}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"FiniteElement", 
        "InterpolationOrder" -> {psiR -> 2, psiI -> 2}, 
        "MeshOptions" -> {"MaxCellMeasure" -> 1}}}][[1]];

Which gives no errors. Then we combine to get the full solution

solR = psiR /. solP[[1]];
solI = psiI /. solP[[2]];
sol[x_,t_] := solR[x,t] + I*solI[x,t];

The problem is, now the "absorbing" boundary condition at xmax is partially reflecting, as can be seen in the animation

Animate[Plot[Abs@sol[x,t], {x, xin, xmax}, 
  PlotRange -> {All, {0, 0.1}}], {t, 0, 800}]

The question is, how should an absorbing boundary condition be properly implemented for this case of coupled PDEs?