# Proper way of applying absorbing boundary conditions to discrete differential equation system

I have 2nd order $$1+1$$ dimensions PDE (a wave equation) that I have rewritten as a coupled system of first order equations and discretized the spatial dimension using a grid of $$N$$ points labeled by an index $$i$$ and spacing $$h$$ such that $$r_i = r_0 + i*h$$ with $$r_0 > 0$$ with the intent to solve it numerically using the method of lines.

To reduce the order of the PDE, I've used the following definitions:

$$\Pi(t,r) \equiv \partial_t\phi(t,r) \quad(1)$$ $$\psi(t,r) \equiv r^2f(r)\partial_r\psi \quad(2).$$

At the endpoint the grid, using fourth order finite difference formulas for the spatial derivatives the system takes the form

$$\partial_t\Pi_N = \frac{f(r_N)}{r_N^2}\left[\frac{25\psi_N - 48\psi_{N-1} + 36\psi_{N-2} - 16\psi_{N-3} + 3\psi_{N-4}}{12h} \right] - \frac{l(l+1)f(r_N)}{r_N^2}\phi_N \quad (3)$$

$$\partial_t\psi_N = r_N^2f(r_N)\left[\frac{25\Pi_N - 48\Pi_{N-1} + 36\Pi_{N-2} - 16\Pi_{N-3} + 3\Pi_{N-4}}{12h} \right] \quad (4)$$

$$\partial_t\phi_N = \Pi_N \quad (5).$$

I would like to apply absorbing boundary conditions at point $$r_N$$, such that at $$r_N$$ the following must be true

$$\partial_t\phi + \partial_r\phi + \frac{\phi}{r} = 0 \quad (6),$$

or equivalently using (1) and (2),

$$\Pi_N + \frac{\psi_N}{r_N^2 f(r_N)} + \frac{\phi_N}{r_n} \quad(7).$$

My question is the following: To apply the abosrbing B. C. to the system, should I:

1. Substitute $$\psi_N$$ from (7) into (3) inside the derivative (the bracketed part of the expression) and substitute $$\Pi$$ from (7) into the derivative part of (4) ?
2. Substitute the whole derivative (the whole bracketed expression) using (6) and then rewrite $$\partial_t\phi$$ using (1) ?
3. Something else entirely?

Anny help would be appreciated.